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ABSTRACT 

An  overview  and  investigation  of  the  more  popular  digital  filter  design  tech- 
niques are  presented,  with  the  intent  of  providing  the  filter  design  engineer  a  com- 
plete and  concise  source  of  information.  Advantages  and  disadvantages  of  the 
various  techniques  are  discussed,  and  extensive  design  examples  used  to  illustrate 
their  application  to  specific  design  problems.  Both  IIR  (Butterworth,  Chebyshev 
and  elliptic),  and  FIR  (Fourier  coefficient  design,  windows  and  frequency  sampling) 
design  methods  are  featured,  as  well  as  the  Optimum  FIR  Filter  Design  Program 
of  Parks  and  McClellan,  and  the  Minimum  p  -  Error  IIR  Filter  Design  Method  of 
Deczky. 
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I.    INTRODUCTION 

Digital  filter  design  today  is  a  rapidly  growing  field  with  myriad  applications 
in  the  areas  of  signal  processing,  image  processing,  filtering,  prediction,  and  estima- 
tion. Extensive  research  has  resulted  in  a  plethora  of  available  information,  which 
can  be  overwhelming  to  the  engineer  presented  with  a  specific  design  problem.  As 
the  title  suggests,  the  purpose  of  this  thesis  is  to  present  an  overview  of  the  more 
popular  design  techniques  in  an  attempt  to  consolidate  the  information  available  in 
the  literature,  and  provide  an  easily  readable  reference  manual.  Detailed  examples 
are  given  to  illustrate  the  application  of  selected  techniques  to  specific  filter  design 
problems. 

Due  to  the  wealth  of  information  available,  a  summary  and  investigation  of 
every  filter  design  technique  is  not  possible.  As  stated  the  more  popular  meth- 
ods are  outlined  in  detail,  however,  references  are  included  as  sources  of  further 
information. 

The  survey  is  threefoM,  consisting  first  of  recursive  IIR  design  methods,  fol- 
lowed by  nonrecursive  FIR  design,  and  finally  computer-aided  design  (CAD).  A 
summary  of  the  chapter  contents  follows  to  enable  the  reader  to  immediately  lo- 
cate the  section  applicable  to  his/her  particular  design  problem. 

Chapter  II,  Recursive  Filter  Design,  presents  Butterworth,  Chebyshev  and 
elliptic  filter  design  from  both  a  traditional  approach,  wherein  analog  prototype 
filters  based  on  design  specifications  are  converted  to  digital  versions  using  the 
bilinear  transformation,  and  a  direct  design  approach  that  eliminates  the  need  to 
determine  an  analog  prototype. 
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Chapter  III,  Nonrecursive  Filter  Design,  includes  an  overview  of  Fourier  co- 
efficient design,  windowing,  and  a  derivation  and  illustration  of  the  method  of 
frequency  sampling.  Among  the  design  examples  presented  is  the  design  of  a  band- 
pass differentiator  and  a  bandpass  integrator. 

In  Chapter  IV,  Computer-Aided  Design,  the  Optimum  FIR  Filter  Design 
Method  of  Parks  and  McClellan  (that  employs  the  Remez  Exchange  Algorithm) 
is  investigated  and  applied  to  a  bandpass  filter  design  problem.  Also  included  is 
a  short  discussion  of  a  method  to  enhance  the  application  of  this  program  to  the 
design  of  high-order  filters. 

To  complete  the  chapter,  the  Minimum  p-Error  Design  Method  for  IIR  Filter 
Design,  that  utilizes  the  Fletcher-Powell  algorithm,  is  presented.  Furthermore,  a 
discussion  of  the  advantages  and  disadvantages  of  these  two  iterative  techniques  is 
included. 
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II.    RECURSIVE  FILTER  DESIGN 

A.  INTRODUCTION 

The  recursive  realization  of  digital  filters  is  advantageous,  in  that  the  desired 
frequency  response  can  be  obtained  using  a  lower  order  filter  than  if  a  nonrecursive 
realization  were  used  (assuming  linear  phase  is  not  required).  This  is  because  the 
filter  frequency  response  is  influenced  by  both  the  poles  and  zeros  of  the  filter 
transfer  function,  whereas  the  frequency  response  of  a  nonrecursive  realization  is 
determined  only  by  the  filter's  zeros. 

Traditionally,  the  design  of  recursive  digital  filters  involves  the  determination 
of  an  analog  prototype  using  one  of  several  methods:  Butterworth,  Chebyshev  or 
elliptic,  to  name  a  few,  and  converting  these  to  a  digital  version  using  impulse- 
invariant  design  or  bilinear  transformation  methods,  the  latter  being  the  more 
popular. 

This  chapter  first  presents  a  review  of  the  traditional  design  techniques  us- 
ing examples  involving  Butterworth,  Chebyshev  and  elliptic  filters;  then  a  Direct 
Design  method  is  presented,  whereby  the  requirement  for  determining  an  analog 
prototype  is  eliminated,  thus  reducing  the  algebraic  complexity  of  recursive  filter 
design. 

B.  TRADITIONAL  DESIGN  TECHNIQUES 

Recursive  filter  design  techniques,  as  stated  in  the  introduction,  involve  deter- 
mining an  analog  prototype  filter  from  given  filter  specifications,  and  then  convert- 
ing the  prototype  to  a  digital  version  using  a  bilinear  transformation.  Since  this  is 
a  well  known  design  procedure  [1],  details  of  the  derivation  will  not  be  considered. 
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However,  as  a  quick  review,  three  examples  involving  the  design  of  Butterworth, 
Chebyshev  and  elliptic  filters  will  be  given. 

It  should  be  noted  that  the  tables  of  analog  prototype  transfer  functions  found 
in  this  section  (Tables  2.2,  2.3  and  2.4),  all  have  a  critical  frequency  of  unity 
(wc  =1).  Thus,  to  obtain  a  filter  transfer  function  with  a  different  critical  fre- 
quency based  on  these  prototypes  requires  the  use  of  an  appropriate  analog  fre- 
quency transformation  (Table  2.1). 

For  example,  suppose  a  lowpass  analog  filter  with  a  critical  frequency  of  wc 
is  desired.  The  following  relationship  applies: 

HLP(jwc)  =  HLPp(jl)  (2.1) 

where  Hlpp  is  the  prototype  filter  with  a  critical  frequency  of  wc  =  1. 

To  convert  the  prototype  transfer  function  to  the  desired  transfer  function, 
Hlp(s),  the  s  in  the  lowpass  prototype  filter  is  replaced  by  s/wc  as: 


HLP(s)  =  HLPp(s] 


(2.2) 


Table  2.1  (Reference  1)  lists  the  appropriate  frequency  transformations  for  high- 
pass,  bandpass  and  bandstop  filters.  A  summary  of  the  design  procedure  follows. 
1.    Bilinear  Transformation  Design 

From  given  analog  filter  specifications,  determine  the  approppriate  s- 
domain  design  technique  (i.e.,  Butterworth,  Chebyshev  or  elliptic). 

Design  the  analog  prototype  filter  in  accordance  with  the  technique  se- 
lected above.  This  involves  the  translation  of  the  filter  specifications  to  those  of 
a  lowpass  prototype.  The  lowpass  prototype  is  then  designed  according  to  the 
translated  specifications. 
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Transform  the  lowpass  prototype  to  the  desired  lowpass,  highpass,  band- 
pass, or  bandstop  filter  according  to  Table  2.1. 

Transform  the  analog  filter  transfer  function  to  a  digital  version  using  the 
bilinear  transformation. 


H(z)  =  H(s) 


(2.3) 


Example  2.1    Butterworth  Filter 

Design  a  digital  filter  that  is  flat  in  the  passband  from  0  to  the  3  dB  cutoff 
frequency,  fc ,  of  2  kHz.  For  frequencies  greater  than  4  kHz,  the  attenuation  should 
be  at  least  10  dB.  The  sampling  r~.tc  is  20  kHz. 

Step  1:      A  Butterworth  design  is  called  for  because  a  flat  passband  is  desired. 
Step  2:      Convert  the  critical  analog  design  frequencies  to  digital. 

Sampling  frequency:      f3  =  20  kHz     =>•     w9  =  40  x  1037r  rad/s 
Cutoff  frequency:  fc  =  2  kHz      ==>     wc  =  4  x  1037r  rad/5 

Stopband  frequency:     fa  =  4  kHz       =>     wa  =  8  x  1037r  rad/5 

9C  =  wcT  =  WJf,  =  {42*™2*  =  0.2*  rad 

Step  3:      Prewarp  the  analog  frequencies  to  yield  the  desired  digital  frequencies. 

w'c  =  tan  (0'c/2)  =  tan(O.lTr)  =  0.325  rad/5 

w'a  =  tan(Q'a/2)  =  tan(0.27r)  =  0.726  rad/5 

Note:      Prewarping  is  done  to  ensure  the  analog  filter  will  ultimately 
yield  a  digital  filter  with  the  correct  critical  frequencies.  This  is  best 
visualized  in  the  following  graph  (Figure  2.1)  of  the  relationship  between 
the  analog  filter  frequencies  and  the  desired  digital  frequencies, 
i.e.,  <=tan(0n/2). 
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TABLE  2.1 
ANALOG  FILTER  FREQUENCY  TRANSFORMATIONS 

(after  Ref.  [1]) 


Filter  Type 


«/-form  s-lorm 

to  find  (to  tind  H{s)  from  tfLP,(s) 

<a\,p  replace  s  in  prototype 

with) 


00 


Lowpass  filter 


Highpass  filter 


Utfjp    —     (Kq 


s2^ 


Butfts 
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Normalizing  to  the  lowpass  prototype,  according  to  Table  2.1: 


w'    =  0.325  — ►  w       =1  rad 
c  c„ 


i'a  =  0.726  — ►  wap  =  0.726/0.325  =  2.234  rad 


CO 


■71/2 


o  e_  e. 


71/2 


Figure  2.1.      Relationship  Between  Analog  and 
Digital  Filter  Frequencies 

Note:  Normalizing  is  done  to  transform  the  prewarped  critical  frequencies 
of  the  analog  filter  to  their  relative  equivalents  in  a  lowpass  prototype  filter  with 
a  critical  frequency  of  wc  =  1.  This  is  done  to  enable  the  designer  to  take 
advantage  of  previously  compiled  tables  or  frequency  response  curves  corresponding 
to  the  transfer  functions  of  these  lowpass  prototypes.  Thus,  a  prototype  filter  can 
be  selected  whose  frequency  response  characteristics  have  the  same  shape  as  the 
filter  that  is  being  designed.  The  transfer  function  for  the  selected  prototype  is 
then  converted  to  a  transfer  function  that  exhibits  the  desired  frequency  response 
characteristics  using  an  appropriate  substitution  (Step  5).  Figure  2.2,  depicts  the 
process  of  normalizing. 
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0) 

0                      „, 
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Figure  2.2.      Normalization 

Step  4:       Find  the  order  N  of  the  lowpass  Butterworth  filter  prototype.  The  order 
N  to  provide  a  gain  of  M^b  at  an  angular  frequency  wa  is 


N 


log10(10(-'vWio)  -1) 

21og10  wa 
log^lO1-!) 


]Md] 


lOdB 


(2.4) 


21og10   2.234 
=  1.37  (Since  N  has  to  be  an  integer  choose,  N  =  2) 
Step  5:      From  Table  2.2  determine  the  transfer  function  for  the  lowpass  filter. 


HLP(s)  =  HLPp(s] 
1 


Hlpp{s) 


s2  +  %/2s  +  1 
0.106 


(.2.5) 


s2  +0.460.5  +  0.106 
Step  6:      Determine  the  transfer  function  for  the  digital  version  of  the  filter  using 

the  bilinear  transformation. 


HLp(z)  =  HLP(s] 


z2  +  2z  +  1 


0.06S(-  +  1)2 


(2.6) 


14.34^  -  U.Oz  +  6.14       z2  -  1.14^  +  0.413 
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TABLE  2.2 
BUTTERWORTH  PROTOTYPE  COEFFICIENTS 

(Table  after  to  Ref.  4) 


N  d\  0,2  <23  a\                    a5                    <^6                   CLj                  <2g 

1  1.0000 

2  1.4141  1.0000 

3  2.0000  2.0000  1.0000 

4  2.6131  3.4142  2.6131  1.0000 

5  3.2361  5.2361  5.2361  3.2361   1.0000 

6  3.8637  7.4641  9.1416  7.4641   3.8637   1.0000 

7  4.4940  10.0978  14.5918  14.5918  10.0978   4.4940  1.0000 

8  5.1258  13.1371  21.8462  25.6884  21.8462  13.1371  5.1258  1.0000 


HLPp(s)  = 


1  +  axs  -f  a2s2  +  . . .  +  aNsN 
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Step  7:      Obtain  a  plot  of  the  digital  filter  frequency  response  to  see  if  the  design 
specifications  have  been  met. 

Figure  2.3  shows  that  the  design  does  indeed  comply  with  the  given  filter 
specifications,  in  that  the  cutoff  frequency  is  0.628  =  0.2-ir  rad  and  the  stopband 
frequency  is  1.257=0.47r  rad  with  gains  of  —3  dB  and  —14  dB,  respectively. 
Example  2.2     Chebyshev  Filter 

Design  a  digital  bandpass  filter  with  the  following  specifications: 

•  1  dB  ripple  in  the  frequency  band  600  to  900  Hz, 

•  sampling  frequency  of  3  kHz, 

•  maximum  gain  of  -40  dB  for  0  <  /  <  200  Hz. 

Step  1:      Convert  the  critical  analog  design  frequencies,  it;,-,  to  the  corresponding 
digital  frequencies,  #,. 

Sampling  frequency:  f9  =  3  kHz 

Lower  ripple  passband  frequency:  fi  =  600  Hz  ==>•  u>£  =  12007T  rad  /  s 
Upper  ripple  passband  frequency:  fu  =  900  Hz  ==>  wu  =  18007T  rad  /  s 
Stopband  frequency:  fa  —  200  Hz  ==>  wa  =  4007T  rad  /s 

0't  =  weT  =  we/fs  =  ^^  =  0.4tt  =  1.26  rad 
0'u  =  wuT  =  wjfs  =  ^^-  =  0.6tt  =  1.88  rad 
B'a  =  waT  =  wa/fs  =  ~~  =  0.1337T  =  1.54  rad 

Ripple  band  center  frequency  :  8'0  =  \/9't9'u  =  0.497T  =  0.418  rad 
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Step  2;      Prewarp  the  analog  frequencies  so  that  the  desired  digital  frequency 

characteristics  will  be  achieved. 

w'(  =  tan(0j/2)  =  tan(1.26/2)  =  0.729  rad/s 

w'u  =  tan (0^/2)  =  tan(1.88/2)  =  1.369  rad/s 

w'a  =  tan (0^/2)  =  tan(0.418/2)  =  0.212  rad/s 


W'Q  =  yjw'ew'u  =  V/(0.729)(1.369)  =  1.00  rad/s 
Ripple  band  is:     B'  =  w'u  -  w't  =  1.369  -  0.729  =  0.64  rad/s 
From  Table  2.1  convert  the  bandpass  filter  frequencies  to  the  correspond- 
ing lowpass  prototype  filter  frequencies. 


wLPp 

= 

B'w'np 

,                                                                (2.7) 

^jp-a.oo)2 

(0.64)w'BP 

WBP            WLPP 

w'e 

0.729           -1 

< 

1.369             1 

< 

0.212         -7.04 

W'Q 

1.000            0 

Step  3:       Determine  the  order  of  the  Chebvshev  prototvpe  filter  that  meets  the 

design  requirements. 

For  a  lowpass  Chebyshev 

filter  the  magnitude  -  squared  characteristic  is 

given  by: 

2 

Hlpp{] 

w) 

—                                                               (n  °} 

-i  +  «»c3K«)                        [-s> 

where  Cyv(w)  is  an  Nth-  order  Chebyshev  polynomial  and  e  indicates  the  degree 
of  ripple. 

In  this  example  the  —40  dB  gain  translates  to  M  =  10-2  for  the  desired 
filter  frequency  response  in  the  stopband,   i.e.,   we  want    llf^'ty)!    <    10-2   or 
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\H(jw)\     <   10~4.    The  design  specification  for  1  dB  ripple  in  the  ripple  band 
corresponds  to  a  value  for  e  of  0.2589. 

Making  an  initial  guess  for  the  filter  order  of  N  =  2  yields: 

2 


H2(jw) 


1 


1  +  0.2589  (2w2a  -1)' 

1 

"  1  +  0.2589  (99.12  -if 

which  is  not  less  than  10-4. 

Proceeding  further  with  N  =  3  yields: 

1 


wa=7.04 

■  =  4.01  x  10" 


(2.9) 


HzUw) 


1  +  0.2589  (4wl-3way 
1 


u>„=7.04 


(2.10) 


2.044  x  10' 


1  + 0.2589(1395.65- 21. 12f 
which  is  less  than  10~4.  Therefore,  the  design  needs  can  be  met  using  a  third-order 
Chebyshev  lowpass  prototype  filter. 
Step  4:     Obtain  the  transfer  function  of  the  third  order  lowpass  prototype  filter 

with  1  dB  ripple.    Looking  at  Table  2.3,  yields  the  following  prototype 

transfer  function: 


Hlpp(s)  = 


0.491 


(2.11) 


s3  +  0.988s2  +  1.2385  +  0.491 

Since  this  is  an  odd  order  filter,  the  constant  0.491  in  the  numerator  was 
selected  to  make  |#(;0)|  =  1.  Recall  that  for  Chebyshev  filter  design,  the  constant 
K  in  the  numerator  of  the  lowpass  prototype  filter  is  selected  on  the  following  basis: 


for  N  odd; 
for  N  even; 


H(j0) 
H(jQ) 


(2.12) 


[1  +  *2] 


211/2 
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TABLE  2.3 
CHEBYSHEV  -  PROTOTYPE  DENOMINATOR  POLYNOMIALS 

(Table  due  to  Reference  4) 
0.5-dB  Ripple  Chebvyshev  Filter  (e  =  0.3493) 


1  *  +  2.863 

2  **  +  1.425*  4-  1.516 

3  *'  4-  1.253**  4-  1.535*  +  0.716  -  (s  +  0.626X5*  4-  0.626*  4-  1.142) 

4  **  +  1.197*3  4-  1-717*2  +  1.025*  +  0.379  -  (,»  +  0.351*  +  1.064X**  4-  0.845*  4-  0.356) 

5  *'  -  1.17251  j*  +  1.9374*'  +  1.3096*2  +  0.7525*  +  0.1789 

(*  4-  0.3623)((*  +  0.1120)2  4.  i.01162][(j  -r  0.2931)2  +  0.62522] 

6  *«  -  1.1592a'  -  2.1718**  -  1.5898*'  +  1.1719*2  +  0.4324*  4-  0.0948 

[(*  -  0.0777)*  -  1.0085*H(*  -  0.2121)*  +  0.7382*][(*  -  0.2898)2  -  0.2702*] 

7  s1  +  1.1512**  4-  2.4126*5  4-  1.8694**  -  1.6479*'  4-  0.7556**  +  0.2821*  +  0.0447 

(*  -  0.2562)[(*  -  0.0570)2  4-  1.0064*][(*  +  0.1597)2  +  0.8001  *][(*  4-  0.2308)*  4- 
0.4479*] 

8  *«  -  1.1461*7  4-  2.6567**  4-  2.1492*5  4-  2.1840**  4-  1.1486*'  4-  0.5736*2  4- 

0.1525*  -  0.0237 

[(*  *  0.0436)*  -t-  1.0050*][(*  4-  0.1242)*  +  0.8520*][(*  4-  0.1859)*  4-  0.5693*] 

[(*  4  0.2193)*  -  0.1999*] 

9  *'  4-  1.1426**  -  2.9027*7  +  2.4293**  -  2.7815*'  +  1.6114**  4-  0.9836*'  + 

0.3408**  -  0.0941*  -  0.0112 

(*  +  0.1984)[(*  -r  0.0345)*  4-  1.0040*][(*  +  0.0992)*  4-  0.8829*]((*  +  0.1520)*  + 
0.6553*][(*  +  0.1864)*  +  0.3487*] 
10  *10  -  1.1401*9  4-  3.1499*»  4-  2.7097*''  4-  3.4409**  -r  2.1442*5  4-  1.5274**  + 

0.6270*'  -  0.2373**  -  0.0493*  -  0.0059 

[(*  -  0  0279)*  -    1.0033*][(*  -  0.0810)*  *  0.9051  *][(*  -  0.1261)*  +  0.7183*] 
[(*  -r  0.1589)*  -  0.4612:][(*  4-  0.1761)*  4-  0.1589*] 
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TABLE  2.3  (Continued) 
CHEBYSHEV  -  PROTOTYPE  DENOMINATOR  POLYNOMIALS 

1.0-dB  Ripple  Chebyshev  Filter  (e  =  0.5089) 

1  *  +  1.308 

2  j:-  0.804*  -r  0.637 

3  *'  -  0.738*1  -  I.02Z*  -  0.327  =  (*  -  0.402H*i  -  0.369*  -  0.886) 

4  **  -  0.716*'  -  1.256*1  -  0.517*  -  0.206  =  is*   -  0.210*  -  0.928)(*i  -  0.506*  -  0.221) 

5  **  -  0.7065**  -  1.4995*'  -  0.6935*1  -  0.4593*  -  0.0817 

(s  -  0.2183)[<*  -  0.0675):  -  0.9735i][(*  -  0.1766)1  -  0.6016-] 

6  *«  -  0.7012*5  -  1.7459**  •  0.S670*'  -  0.7715*2  -  0.2103*  -  0.0514 

[(s   -  0.0470)=  -  0.98I7:][(*  -  0.1283)-  -  0.71S7i]{(*  -  0.1753)i  -  0.2630:] 

7  j7  -  0.6979*'  -  1.9935*'  -  1.0392**  -  1.1444*'  -  0.3S25*:  -  0.16661a  -  0.0204 

(s  -  0.1553)[(*  -  0.0346)^  -  0.9867i][(*  -  0.0968):  -  0.7912:][(*  -  0.1399)-  - 
0.43911] 

8  *•  -  0.6961*7  -  2.2423**  -  1.2117*'  -  1.5796**  -  0.5982*'  -  0.3587*:  - 

0.0729*  -  0.0129 

[(*  -  0.0265):  -  0.9898:][(*  -  0.0754):  -  0.8391  :][(*  -  0.1129):  -  0.5607:] 

[(*  -  0.1332):  -  0.1969:] 

9  *»  -  0.6947*'  *  2.4913*"  -  1.3837**  -  2.0767*'  ~  0.8569**  -  0.6445* J  - 

0.16844*:  -  0.0544*  -  0.0051 

(*  -  0.1206)((*  -  0.0209):  -  0.99|9:][(*  _  0.0603):  -  0.8723:][(*  -  0.0924): 
-r  0.6474:][(*  -  0.1134*:  -  0.3445:) 
10     *'»  -  0.6937*'  -  2.7406*»  -  I.5557*7  -  2.6363**  -  1.1585*'  -  1.0389**  + 
0.317S*'  -  0.1440*:  -  0.0233*  -  0.0032 

[(*  -  0.0170)1  -  0.9935:][(*  -  0.0767):  _  0.7I13:][(*  ^  0.0493):  -  0.8962:] 
[(*  -  0.0967):  _  0.4567:][f*  -  0.1072):  -  0.15741] 


2-dB  Ripple  Chebyshev  Filter  (e  =  0.7648) 

1  *  +  1.965 

2  *i  -  1.098*  -  1.103 

3  *'  -  0.988*1  -  1.238*  -  0.491  =  (*  ->-  0.494H*2  -  0.490*  -  0.994) 

4  **  -  0.953*'  -r  1.454*1  -r  0.743*  ■*■  0.276  =  (*i  -  0.279*  -  0.987)(*i  -  0.674*  -  0.279) 

5  *'  -  0.9368**  -  1.6688*'  -  0.9744*1  -  0.5805*  4-  0.122? 

(*  -  0.2895)[(*  -  0.0895)1  ■*-  0.99O1  ■}[{:  *  0.2342)1  -  0.61191] 

6  **  a.  0.9282*'  -r  1.9308**  -  1.2021*'  -  0.9393*1  -u  0.3071*  -  0.0689 

[(*  -  0.0622)1  -  0.9934i][(*  -,-  0.1699)1  j.  0.7272=][(j  -  0.2321  )i  -  0.26621] 
*7  -  0.9231*'  -  2.1761*'  -  1.4288**  -  1.3575*'  -  0.5486*1  -  0.2137*  -  0.0307 
(*  -  0.2054)((*  -  0.0457):  -  0.9953i]K*  -  0.1281)1  -  0.7982i][(*  -  0.1851 )-  - 
0.44291] 

8  *•  -  0.9198*7  -  2.4230**  -  1.6552*'  -f  1.8369**  -  0.8468*'  -  0.4478*1  + 

0.1073*  -  0.0172 

[(*  -  0.0350)1  -  0.9965i][(*  -  0.0997)1  -  0.8448i][(*  J.  0.1492)1  -  0.56441] 

[(*  -  0.1759)1  -  0.19821] 

9  *»  -  0.9175*«  -  2.6709*'  -  1.8815**  -  2.3781*'  -  1.2016**  -  0.7863*'  4. 

0.2442*1  -  0.0706*  -  0.0077 

is  -  0.1593)[(*  -  0.0277)1  _  o.9972i][(*  -  0.0797)1  +  0.8769i][(*  -  0.1 221  )•  - 
0.6509i][(*  -  0.1497)1  -  0.34631] 
0  j'o  -  0.9159*»  -  2.9195*»  -  2.1079*7   -  2.9815**  -  1.6830*'  -  1.2445**  - 

0.4554*'  -  0.1825*1  -  0.0345*  -  0.0043 

[(s  -  0.0224)1  -  0.9978i][f*  -  0.1013)1  -  0.7143i)[(*  ~  0.0651)1  -  0.9001 :] 
[(*  -  0.1277)1  -  0.4586!)[(*  -  0.1415)1  -  0.15801] 
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Step  5:     Convert  the  lowpass  filter  prototype  to  a  bandpass  filter  prototype. 
From  Table  2.1: 


HBp{s)  =  HLPp{s) 


(2.13) 


(oft)    +  0-9*8  (£±I)    +  1.238  (&£)+ 0.491 

0.129s3 

~  s6  +  0.632s5  +  3.507s4  +  1.394s3  4-  3.507s2  +  0.6325  +  1 

Step  6:     Determine  the  digital  version  of  the  transfer  function  using  the  bilinear 

transformation. 


HBP(z)  =  HBP(s) 


-n3 


0-129  (frf) 


({=£)•  +  0.632  (f-i)5  +  3.507  (-1)4  +  1.394  (f-1)3  +  3.507  ({£)' 

-(-  1.238  I  ^"— T  )  +0.491 


x2 

0.011  (z6  -  3z4  +  3;2  -1) 


z6  +  2.15324  +  1.786z2  +  0.545 

(2.14) 

Step  7:  Confirm  that  the  filter  design  meets  specifications  by  obtaining  a  fre- 
quency response  plot.  Figure  2.4  indicates  the  specifications  for  a  pass- 
band  of  1.26  to  1.88  rad  with  1  dB  ripple  have  been  met. 

Example  2.3      Elliptic  Filter 

A  lowpass  elliptic  filter  with  the  following  specifications  is  desired: 

•  passband  ripple  of  0.5  dB, 

•  passband  ripple-edge  frequency  of  2  kHz, 

•  stopband  gain  should  be  at  most  —20  dB  for  frequencies  greater  than  4  kHz, 
and 

•  sampling  frequency  is  20  kHz. 
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Step  1:     Convert  the  critical  analog  design  frequencies  to  digital  frequencies. 

Sampling  frequency:  f3  =  20  kHz 

Passband  ripple-edge  frequency:      /i  =  2  kHz    ==>   w \  =  4  x  1037r  rad  / s 

Stopband  ripple-edge  frequency:     f2  =  4  kHz    =>  w2  =  8  x  1037r  rad  Js 

(4  X   103)  7T 

0[  =  WlT  =  wylfa  =  V20xlQ;     =  0.2;r  =  0.628  rad 

(8  X  103)  7T 

9'2  =  w2T  =  w2/f3  =  V20xlQ3     =  0.6tt  =  1.88  rad 

Step  2:  Prewarp  the  analog  frequencies  in  order  to  determine  the  appropriate 
lowpass  prototype  filter. 

w[  =  tan(0j/2)  =  0.325  rad/s 
w'2  =  tan  (02/2)  =  1.376  rad/5 

Step  3:  Determine  the  order  of  the  elliptic  prototype  filter  that  meets  the  design 
specifications.  Find  the  ratio,  R,  of  the  stopband  frequency,  w'2,  and  the 
passband  frequency,  w[ . 

R  =  u»iM  =  J||=  4.234  (2.15) 

From  Table  2.4b  determine  the  prototype  filter  order  and  the  resulting 

transfer  function.    Looking  at  Table  2.4,  it  can  be  seen  that  a  filter  of 

order  N  =  2  has  a  value  for  R  of  2.76261,  which  is  more  than  sufficient 

to  meet  the  design  specifications. 

Step  4:     Obtain  the  transfer  function  of  the  second-order  lowpass  prototype  filter. 

Again,  looking  at  Table  2.4,  the  prototype  transfer  function  is  determined 

to  be: 

_       ,  ,  0.1s2  +  0.5338  ,__, 

Hlp'{3)  =  ,2+0.80943  +  0.5667  (2-16) 
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TABLE  2.4a 

ELLIPTIC  FILTERS 

Generalized  Transfer  Functions 

(due  to  reference  5) 


iV  =  2 


-V 


Ho  s2+A0l 


s  +  s0     s2  +  Bns  +  Boi 

H0S2  +H0A01 


s3  +  (BU  +s0)s2  +(B0i  +Buso)s  +  Bolso 


iV  =  4 


(s2  +A01)  (s2+A02) 

s2  +Bus  +  B01  '  s2  +  Bus  +  B02 

H0s*  +Hq(A01  +A02)s2  +HqAq1Ao2 

s4  +  (5„  +  B12)s3  +  (B02  +  BnB12  +  B0i)s2  +  (5nB02  +  B0iB12)s 

+  Bq\Bo2 


H4(s)  =  H0 


H4(s) 


N  =  5 

#0s4  +  #o(.4oi  +Ao2)52F0AoiA( 


#*(*)  = 


s5  +(5„  +  £12  +s0)64  +  [BQ2  +  BUB12  +B01  +Buso+B12s0]s* 
+  [BUB02  +B0iB12  +B02s0  +BuB12so  +  B01so]s2 


+  [BoiBq2  +  BhBo2Sq  +  I?oi#12<So]  5  +  BoiB02Sq 
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TABLE  2.4b 
ELLIPTIC  FILTERS 

(due  to  reference  5) 


n- 


Passoand  ripple  0  5  dB.  stopoand  gains    -20.  -30.   -40.   -50. 
Passoand  nppie  «  0  5  dB.  stopoand  gain  • 


■60.  -70  i 
-20d8 


33789 

"5640 


6507ft 
0721 1 


0  566660 
0  808321 


i)  8277X7 
0  973640 

0  61 1899 

0  934X30 
0.990620 


0  97h474 
0  94668  I 


00775 
00 1  OX 


0 809390 
I)  359160 


0  931959 
0  1 3A543 

0.41281ft 

0  049395 

0.933855 
0.156221 

0017576 

0  413652 
0  0563X4 
0  006219 

0  933X64 
'I  156548 

0  020051 
0  002197 

0  4  I <695 
0  056505 
0.OO7O93 

0  000775 


306214F. 

667292 


303X95E  • 
667292 


I  303X6 IF. 
I  667292 


0  I00I92E  +  OO0 


0  3037X6E-OO0 


2.76261 
1.421X9 

I.I3IX8 

1.04465 

I  01553 

1.00545 


Passoand  rippie  =  0  5  dB;  stopoand  gain 


0  318702 
0  5973X4 


0  639007 
0  3X2044 


I6294E 
21X7XE- 


808X0 
.92322 


3X6X0 
07474 

!  13439 
15171 

03I6X 
,  <794l 
3X394 
06301 
01359 

:  13409 
1 5070 


0  398996 
0  798764 

0  64X724 
0  907216 
0  402050 


0.650591 
0  9207X5 
0  981925 


0  964X98 
0  992137 

0 650636 
0  921256 
0  984652 
0  996<X9 


0  822201 
0  191032 
0  480774 
0  0X80X0 
0  828822 
0  237025 
0.039181 

0  484325 
0  108419 
0.017159 


0  04804  3 
0  007464 

0  484454 
0  109144 
0  021005 
0  003238 


XX07E  - 
1761 
6296E  - 


0  ||8701E*000 

0  511761 

0  3I62S9E-OOI 


12912 
05394 
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TABLE  2.4c 
ELLIPTIC  FILTERS 

(due  to  reference  5) 


PassDand 

rippie   1  dB 

.  stopDand  gains 

-20.  -30. 

-40.  -50.  -60.  - 

70  dB 

PassDand  ripple 

=  1  dB:  siopband  gain  =   -  20  dB 

N 

Aa 

Ba 

S„ 

Ho-S,, 

ft 

: 

i 

4  4  2.14  J 

0  497233 

0  676727 

0  1001851-  -000 

2  32474 

» 

1 

1  5HS65 

0.790229 

0.282927 

0  2X1  (WOE-*  000 

0  565168 

1   30797 

* 

3  81475 
1    159S6 

0536633 
0926578 

0  768217 
0  099029 

0  IOOI85L  +  000 

1  09029 

' 

;. 

1  5IH52 

1    IU886 

0  K0XO49 

0  975703 

0  318242 
0  0.12771 

0  279829b  -000 
0  565168 

1  02826 

- 

i 

.1  80873 

1   MIA? 
1  015V) 

0  537071 
0  931 194 
0992101 

0 769217 

0  1 10761 
0.010654 

0  I0OIX4K-0O0 

1.00902 

7 

; 

1  51782 
1  04421 
1  (XU97 

0.808242 
0.977941 
0  997446 

0.318627 
0  036573 
0.003444 

0  2798I4E-0O0 
0.565168 

1.0O29O 

" 

5 

3  80866 
1  14350 

1   014115 
1  00160 

0.537076 
0  933266 
0  992834 
0  999176 

0  769228 
0  110888 
0.01 1881 
0.001  III 

0.I0OI85E-OO0 

1  CXXN3 

' 

i 

4 

1  51X12 
1  04421 
1  IKM5I 
1  00052 

0  808099 
0  977936 
0  997680 
0  99*1734 

0  118714 
0  036644 
0  003845 
0  000.159 

0.27961  It  -000 
0  565 16* 

1  000.30 

Passband  ripple 

*   1  dB,  stppband  gam  -    -  30  dB 

N 

Aa 

So 

8„ 

H0iS0 

11, 

2  2029.1 

0  586596 

0  31 1981 

0  1 1  3223E 
0  430700 

+  000 

1  73254 

5  72672 

0  348  1  1  1 

0  682880 

03162861": 

-001 

1  25040 

1  37628 

0  803477 

0  148347 

1  96914 

0  63.1529 

0  383990 

0  1 1 1 1 80E 

toon 

1  09554 

1    11918 

0  914405 

0  0646 12 

0  430700 

5  67743 

0  350184 

0  687075 

0  3I6278E 

•001 

1  0.1799 

1   11626 

0  828467 

0.179733 

1  05468 

0  964084 

0  0271  II 

1  96122 

0  614884 

0  186055 

0  III  1231. 

-IKK! 

1  01536 

1    1 1 870 

0  925872 

0  077672 

0  430700 

1  02200 

0  985162 

0  01  1201 

5  67614 

0350237 

0  68  7189 

0  3I6265R 

-001 

1  00625 

1   31463 

0829171 

0  180623 

1  04690 

0  969003 

0032478 

1  00893 

0  993908 

0  0O459X 

1  96311 

0  614904 

(I  386117 

0  IIII09E- 

.  (K)0 

1  00255 

1    1 1814 

0  926187 
0  987211 

0  078043 
0  013399 

0  430700 

1    <X)364 

0  997505 

0  001883 
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TABLE  2.4d 
ELLIPTIC  FILTERS 

(due  to  reference  5) 


Pwacano 

nople;  2  a8 

slopoand  gams 

-20.   -30. 

-40.   -50.     -60.    - 

70  dB 

Passoana  ripple 

»  2  dB:  sioooand  gam  =   -  20  dB 

N         i 

Aa 

So 

a„ 

H0/S0 

11. 

:         i 

<  60961 

0  454891 

0  5.17126 

0  1001036 -MXX) 

1  94332 

'         ' 

1   429.19 

0  74 .11X0 

i)  204089 

0.25444 IE  4-000 

II  458848 

1  20808 

*    ! 

1   10765 

0  486218 
0  4)5564 

U  597266 
11  06.1585 

O.IOUI02E  +  000 

1  05569 

5    : 

1   191  16 
1  02976 

0  807116 
0  481070 

ii  22.1995 
I)  018680 

0  2538786 +  IXX) 
0  458848 

1  01567 

«    \ 

i  25657 

1  IIWI  » 

i  i«)845 

0  486417 
II  940286 

0  44JS12 

I)  597679 

1  1)694  17 
1)  005.14ft 

0  IOOI02E  +  000 

1  (XW47 

1   "0242 

0  807184 

ii  020.162 

o  25.1X.V)E*(XX) 
II  458848 

1  00128 

j 

'  2S73A 

1  0441  1 

i  (x)7s: 

!  1)0069 

0  486114 
0  940275 
0,994418 
0  994548 

0 597662 
0  0694X8 
0. 00588.1 
0  000446 

1)  1000406*000 

1 .00037 

; 

l  IP7<5 

0  807065 

o  224111 
0  020419 

0  25.141X6*000 

1  0001 1 

Passoand  riooie   =  2  dB.  siopoand  gain  =    -  30  dB 


41395.1 

0  3I6259E 

238474 

0  1041836 
0  345928 

517  241 
105666 

o  .11625X6 

2XS64I 
042702 

i)  I02947E 

514562 
124605 
016632 

0  3162596 

286656 
050009 

o  1029206 
II  345928 

050171 
007452 

00119  10 
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Step  5:  Find  the  analog  transfer  function  of  the  desired  filter  based  on  the  lowpass 
prototype. 

This  is  done  by  finding  a  scaling  factor,  a,  to  convert  the  lowpass  proto- 
type's cutoff  frequency,  wip,  to  the  desired  cutoff  frequency,  w\. 

For  the  prototype: 

fl  =  2.76261     ;     ^  =  1/^  =  0.6016 

For  the  desired  filter: 

wl  =  0.325 


Therefore,  the  scaling  factor,  a,  is: 

0.6016 


a  =  wip/wi 


0.325 


1.85 


(2.17) 


Thus,  the  actual  transfer  function,  Hlp(s),  can  be  found  as  follows: 


HLP(s)  =  HLPp(s] 


Rlpp{*)\ 


0.1s2  +  0.5338 


(2.19) 


s2  +  0.80945  +  0.5667 
_  0.1(1.85s)2+ 0.5338 

~  (1.85s)2  +0.8094(1.855)  +  0.5667 

0.3423s2  +  0.5338 
"  3.423s2  +  1.497s +  0.5667 
Step  6:    Determine  the  digital  version  of  the  transfer  function  using  bilinear  trans- 
formation. 


Hlp(z)  =  HLP(s) 


0.3423(2-1/2  +  1)'  +0.5338 


3.423(2  -  1/2  +  iy  +  1.497(2  -  l/x  +  1)  +  0.5667  (2.20) 
0.876122  +  0.3832  +  0.8761 


5.48722  -  5.71262  +  2.493 

0.159722+ 0.06982 +  0.1597 

22  -1.04112  +  0.4543 
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Step  7:    Verify  that  the  design  meets  the  specifications  from  a  frequency  response 
plot  (Figure  2.5). 
2.    Summary 

This  concludes  the  review  of  the  use  of  the  bilinear  transformation  to 
design  recursive  filters  based  on  Butterworth,  Chebyshev  or  elliptic  analog  pro- 
totypes. As  can  be  seen,  this  method  is  very  involved  in  terms  of  algebraic  ma- 
nipulation. The  following  sections  will  introduce  the  direct  design  technique.  It 
will  be  shown  that  the  direct  design  procedure  reduces  considerably  the  amount  of 
algebraic  calculations  required,  and  eliminates  the  somewhat  confusing  procedure 
of  prewarping. 

C.    DIRECT  DESIGN  TECHNIQUE 

The  direct  design  technique  is  based  on  the  bilinear  transformation  in  the 
following  manner: 

The  tables  of  coefficients  for  analog  lowpass  prototype  filters,  i.e.,  Butterworth  - 
Table  2.2,  Chebyshev  -  Table  2.3,  and  elliptic  -Table  2.4,  were  converted  to  z- 
domain  versions  through  the  use  of  the  bilinear  transformation,  and  can  be  shown 
to  have  a  critical  frequency  of  7r/2.  The  reason,  as  stated  in  Section  A,  is  that 
the  analog  prototype  filters  comprising  these  tables  all  have  a  critical  frequency  of 
wc  =  1.  For  the  bilinear  transformation,  the  following  substitution  for  s  is  used  in 
the  analog  prototype  transfer  function  to  convert  it  to  a  digital  prototype  transfer 
function: 

-St  «  • =B  w 

An  analog  critical  frequency  of  5  =  jwc  =  jl  implies: 

2=1-±l^  =  l±l  (2.22) 

l-jwc        I- j 
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c 

—  c 

=5  d 


~^--l           !           !    •       i 

\          ! 

I  \        ! 

I    \ 

i 

\                 ! 

\    1            ! 

\i 

i      x              i 

j             j\          !             ! 

k                 ^  .l- — ~~~T~ 

0.00000   0.62832      1.25G64       J. 88-196 

FREQUENCY,  THETA 


2.51328 


3.14160 


N  = 


Figure  2.5.     Frequency  Response  for  Elliptic 
Lowpass  Filter 
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since  z  —  ejB ,  this  gives: 

e*°c  =  1±4  =  leJ7r/2  (2.23) 

Thus,  it  can  be  seen  that  an  analog  critical  frequency  of  wc   =   1  yields  a 

digital  critical  frequency  of  6C  =  ir/2.  Tables  of  prototype  filters  were  thus  created, 

specifically,  Butterworth  -  Table  2.7,  Chebyshev  -  Table  2.8  and  elliptic  -  Table  2.9. 
The  use  of  these  tables  to  design  a  digital  filter  based  on  analog  specifications 

involves  the  following  generalized  steps.     To  further  illustrate  the  technique,   a 

detailed  description  of  its  derivation  and  use,  including  examples  for  the  various 

filter  types,  will  be  presented. 
1.    Direct  Design 

Step  1:  Determine  the  appropriate  desired  filter  type  based  on  the  design  speci- 
fications, i.e.,  Butterworth,  Chebyshev  or  elliptic. 

Step  2:     Convert  the  analog  filter  frequency  specifications  to  digital  equivalents. 

Step  3:  Use  the  digital-digital  frequency  transformations  of  Table  2.5  to  normal- 
ize the  desired  filter's  design  frequencies,  so  that  the  appropriate  lowpass 
prototype  filter  may  be  selected. 

Step  4:    From  the  normalized  design  frequencies  obtained  in  Step  3,  determine 
the  lowpass  prototype  that  meets  or  exceeds  these  requirements. 
For  Butterworth  and  Chebyshev  filters  -  this  is  done  using  the  design 
curves  of  Figures  2.8  and  Figures  2.10-2.12. 
For  elliptic  filters  -  Table  2.9  is  used. 

Step  5:  Obtain  the  actual  filter  transfer  function  from  the  lowpass  prototype, 
based  on  the  digital-digital  frequency  transformations  of  Table  2.5. 
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TABLE  2.5 
DIGITAL  FILTER  FREQUENCY  TRANSFORMATIONS 

(Due  to  Reference  2) 


Type  Transformation 

(Replace  z  in  LP  digital 
prototype  with) 


Design  Constants 


1. 

Lowpass 

i^t 

_  sln(ec/2-0'c/2) 
sin(9c/2+9'c/2) 

2. 

Highpass 

_   2-a 
1  —  az 

cos(9c/2-9'c/2) 
01         cos(9c/2+9'c/2) 

3. 

Bandpass 

k: 

cos(0'u/2+^/2) 
01  ~  cos(9'u/2-9't/2) 

=  tan  ^cot(0'u/2- 0^/2) 

4. 

Bandstop 

z2  - -22- z+ ±=Jl 

k  = 

cos(9'u/2+9'(/2) 
~  cos(9'u/2-9'(/2) 

=  tan^tan(0'u/2-0'€/2) 
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As  can  be  seen  from  the  previous  design  summary,  the  crux  of  this  method 
is  the  use  of  the  digital- digital  frequency  transformations  of  Table  2.5,  which  merit 
explanation  [2]. 

These  transformations  enable  the  user  to  convert  the  lowpass  digital  pro- 
totype, HlPp(z)  (be  it  Butterworth,  Chebyshev  or  elliptic)  to  the  actual  lowpass 
(LP),  highpass  (HP),  bandpass  (BP),  or  bandstop  (BS)  filter  transfer  function. 

The  transformations  provide  a  means  of  transferring  the  stability,  inherent 
in  the  prototype  filter,  to  the  actual  filter;  that  is,  the  poles  of  the  actual  filter  will 
lie  inside  the  unit  circle,  as  they  do  for  the  prototype  filter.  For  this  reason  the 
frequency  response  of  the  lowpass  prototype  filter  at  a  specific  frequency  of  0a  must 
be  the  same  as  the  desired  value  of  the  frequency  response  of  the  actual  filter  at 
its  corresponding  frequency  of  6'a . 

HLPp(e^)=HType(e^)  (2.24) 

where, 

0a  =   prototype  frequency 

9'a  =   desired  filter  frequency 

Type  =  LP,  HP,  BP,  or  BS 
For  example,  when  transforming  a  lowpass  prototype  filter  to  a  desired 

lowpass  filter,  one  wishes  to  maintain  the  integrity  of  the  magnitude  characteristic 

of  the  prototype  filter,  while  expanding  or  compressing  the  frequency  scale  so  that 

the  critical  frequency  is  changed  from  the  prototype  value  of  8C   =  7r/2  to  the 

critical  frequency  of  the  actual  filter  B'c,  (see  Figure  2.6). 

According  to  Table  2.5,  this  involves  the  following  substitution  for  z  in 

the  lowpass  prototype  transfer  function. 

z  =  ±=2-  (2.25) 

1—0:2; 
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^ 

Lowpass 

Prototype 

<^M 

Lowpass 

1.0 

-^^^ 

1.0 

\^ 

0.707 

\ 

\ 

,\ 

v_ 

0 

8e 

K 

6 

9 

e 

Figure  2.6.      Lowpass  Prototype  -  Lowpass  Transformation 

For  the  critical  frequency  of  the  prototype,  9C,  to  map  to  the  critical 
frequency  of  the  actual  filter  &c ,  the  following  must  be  true: 

Substituting,  e*9c  —  z  on  the  left  in  Equation  (2.25),  and  eJ  c  =  z  on  the 
right  gives, 


e>°c  = 


eJI 


1  -  aeJ"c 
and  solving  for  the  design  constant,  a,  yields, 

sm(9c/2-0'r/2] 


(2.26) 


sin(*c/2  +  *c/2) 

Similar  arguments  apply  to  HP,  BP,  and  DS  filters. 

Thus,  the  transformations  of  Table  2.5  accomplish  the  appropriate  map- 
ping from  the  lowpass  prototype  digital  fiiter  to  the  actual  filter,  and  ensure  that 
stabilty  is  maintained. 
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2.    Butterworth  Filters  -  A  Direct  Design  Approach 

The  design  curves  shown  in  Figures  2.8a  -  d,  describe  lowpass  prototype 
filters  with  a  critical  frequency  of  9C  =  7r/2,  and  of  order  N  =  1  through  N  =  8. 
Higher  order  filters  can  be  realized  by  cascading  the  appropriate  lower  order  filters 
characterized  by  these  curves. 

The  curves  were  constructed  by  using  the  s-domain  version  of  the  But- 
terworth filter  coefficients  for  N  =  1  through  N  =  8,  and  then  determining  the 
2-domain  version  of  these  transfer  functions  through  the  use  of  the  bilinear  trans- 
formation and  the  conversion  tables  obtained  from  Reference  3,  (see  Table  2.6). 
These  tables  establish  the  necessary  relationships  between  the  coefficients  of  an 
analog  filter  and  its  digital  counterpart. 

a.    Generalized  Analog  Transfer  Function 


AQ  +Als  +  A2s2  +  ...  +  Aksk 
BQ+B1s  +  B2s2  +...  +  Bksl 


h(s)  =  ?  t„a:t  ?:,:•••;„'%  (2.27) 


b.    Generalized  Digital  Transfer  Function 

m  s       ap+aiz-1  +a2z-2  +...  +  akz-k 

H{Z)  ~    l  +  61*-i+62z-»+... +  ***-*  (2'28) 

For  example,  in  the  case  of  N  =  1  the  bilinear  transformation  digital 
filter  coefficients,  in  terms  of  the  analog  filter  coefficients  are: 

A  Bo+Bi 

a0  (A0+Ai)/A 

ai  (Ao-Ai)/A 

h  (Bo-B^/A 

Note:  C  is  the  critical  frequency,  wc  =  1. 

For  a  first-order  Butterworth  prototype  filter  described  by: 

H(s)  =  -L-  (2.29) 

5  +  1 
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TABLE  2.6 

BILINEAR  TRANSFORMATION  DIGITAL  FILTER  COEFFICIENTS 

IN  TERMS  OF  ANALOG  FILTER  COEFFICIENTS 

(due  to  reference  3) 


A      Bo+8,  OS^C^ejC3 


Bo*B,C 
(AQ*A,C)/A 
(Aq-A^J/A 
(80-e,C)/A 


(1st  order) 


80*8,0*820^ 
(a0*a,c*a2c2)/a 

(2A0-2A2C2)/A 
(A0-A,C*A2C2)/A 
(280-2B2C2l/A 
(8o-8,C*82C2)/A 

(2nd  order). 


(A0«-A,C*A2CZ+A,Ci)/A 


Bo-8iC*32C2~8jC3«84C4 

(A0-A.C*A2C2-A3C3'-A<,C4)/A 
(4A0*2A, C-2A jC3 -«A„C4 ) /A 
(6A0-2A2C2»6A4C4)/A 
(4A0-2A,  C-2AjC3-»A4C4  :  /A 

(A0-6,c+A2c2-a3C3-Aac4)/A 
(1B0*2B,C-293C3-484C4VA 
(680-2B2C2*684C4)/A 
<4B0-2B,C*2B3C3-«J8,.C4)  /A 
(8o-a,C*82C2-&3C3*e4C*}/A 

(4th  order) 


( 3A0*A, C-A2C2-3A jC3 ) /A 

(3A0-A,C-A2C2*3A3C3)/A 

<A0-A,OA2C2-A3C3)/A 

(380*8,C-B2C2-3BjC3)/A 

(SBq-B,  C-82C2i-383C3)/A 

(B0-81C+B2C2-B3C3)/A 

(3rd  order). 

80*81C*B2C2*ejC3*e4C4*B,C9 

( AQ*A, C* A2C2* A3C3 *A4C4*A,C5 ) /A 

(5A0*3A.C+A2C2-A3C3-3A4C4-5A5C5)/A 

{ I0A0+2A, C-2A2C2-2A3C3  *2A4C4 ♦ lOA,C5 ) /A 

(IOA0-2A,C-2A2C2*2A3C3+2A4C4-IOA5C5!/A 

(5A0-3A1C*42C2*A3C3-3A4C4*5A5CS)/A 

(A0-A,C*A2C2-A3C3*A4C4-A5C5)  /A 

(5Bo*38,C+e2C2-BjC3-3B4C4-5S5C5)/A 

( 1003*28, C-2e2C2-28jC3*2B4C4*iOB5C5) /A 

(iC80-28,C-282C2*28jC3*284C4-iOe5C5UA 

(580-3B,  0*82^*8 jC3-3B4C4*5B5C,)/A 

(8o-8,C*82C2-83C3*B4C4-B5C5)/A 

(5th  order) 


the  values  for  the  digital  filter  coefficients  are: 

A  1+1=2 

a0  (l+0)/2  =  0.5 
aj  (l-0)/2  =  0.5 
61  (1  -  l)/2  =  0 

The  transfer  function  H (z)  is  therefore: 

0.5  +  0.5Z"1 


Hi(z) 


Hi(z) 


0.5z  +  0.5 


(2.30) 


1  z 

A  plot  of  this  transfer  function  is  illustrated  by  the  curve  (N  =  1)  in  Figure  2.7 
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Table  2.7  is  a  summary  of  the  analog  and  digital  versions  of  lowpass  prototype 
Butterworth  filter  transfer  functions  for  filter  orders  up  to  N  =  8. 

To  illustrate  the  direct  design  approach  the  design  problem  of  the 
previous  section  will  be  used  (Example  2.1),  where  the  problem  statement  was: 

Design  a  digital  filter  for  a  20  kHz  sampling  rate  that  is  flat  in  the 
passband  of  0  to  the  3  dB  cutoff  frequency  of  2  kHz,  and  has  a  gain  of  not  more 
than  —10  dB  for  frequencies  greater  than  4  kHz. 

Step  1:     A  Butterworth  design  is  called  for  because  a  flat  passband  is  desired. 
Step  2:     Convert  the  critical  analog  design  frequencies  to  digital. 

Sampling  frequency:      fa  =  20  kHz 

Cutoff  frequency:  fc  =  2  kHz      =»     w c  =  4  x  1037r  rad/s 

Stopband  frequency:     fa  =  4  kHz       =*►     wa  =  8  x  1037r  rad/s 

ffc  =  wcT  =  »„//.  =  {42*™*l*  =  0.2*  rad 
ffa  =  WaT  =  wa/f.  =     20XX101Q3     =  0.4t  rad 

Step  3:  Since  the  "design  chart"  is  based  on  a  critical  frequency  of  9C  =  t/2, 
digitized  design  specifications  need  to  be  normalized  to  this  frequency, 
so  that  the  order  of  the  required  lowpass  prototype  filter  can  be  deter- 
mined. For  this  process  the  following  frequency  transformations  apply 
(from  Table  2.5,  prime  denotes  desired  filter). 
•  Prototype  filter  from  desired  filter: 


eJ*»  —  a 

~, Tor  (2-31) 

1  -  oe^«  v         ' 


44 


TABLE  2.7 

LOWPASS  PROTOTYPE  BUTTERWORTH  FILTER 

TRANSFER  FUNCTIONS 

(due  to  reference  4) 


Filter  Order 

K  HLPp(s)  HLPp(z) 


1 

3+1 

1 

0.5(2+1) 

z 

(z  +  D2 

a2  +  1.41423+l 

1 

3.41z2+0.59 
(z  +  1)3 

s3+2s2+2s+l 

1 

6z3  +  2z 

(z+D* 

34+2.613l33  +  3.4142s2+2.6131 
1 

9+1 

10 

.6383z4+5.17z2 
(z+1)5 

35+3.236l34+5.2361a3  18. 94z5  +  l  1.996z3+ 1.05492 
+5. 2361a2+3. 23613+1 

i ii+nl 

s6  +3. 8637a5  +7. 4641a4  +9. 141  6a3  +7. 4641a2  33  .7972 6  +2 6.284  z4  +3  .86  22+0  .05  92 

+3.86373+1 

1 (*  +  l)7 

a7+4.4940  36  +  10.097835+14.5918a4  60.367z7  +  55.537z5  +  11 .63323  +  0.46  24  2 

+  14.591833  +  10.0978  32+4.4940s+l 

1 (*  +  l)8 

a8+5.1258a7  +  13.1317a6+21.846235  107.89628  +  1 14.438z6+31 .496  z4 

+  25.688434+21.846233+13.137l32  +2.162z2+0.008 
+5.12583+1 
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•  Desired  filter  from  prototype  filter: 

ejd'«  +  a 


1  +  aeie* 
where  (2.32) 

_sin(flc/2-fl'c/2) 
01      sin(0c/2  +  0'c/2) 

0'c   and  dj,  are  the  critical  and  stopband  frequencies  (respectively)  used  to 
determine  the  optimum  prototype  filter  order. 

Equation  (2.31)  is  applicable  to  this  problem,  and  thus: 

0a  =? 

0'a   =  0.47T 

9C  =  0.5tt 

d'c  =  0.2;r 

_  sin(0.57r/2-0.27r/2) 
a~  sin(0.57r/2  +  0.27r/2) 
=  sin(0.157r)  =  0.454  = 
sin(0.357r)       0.891 


Therefore, 


>j0a  _    eJ°Aw  -  0.510 
"  1  -  0.519e>0-4ff 
=  0.97e^779   =    .2.349 
0.97e->057 

=>         6a  =  2.349  rad   =  0.75tt  rad 


(2.33) 


Step  4:    Once  the  normalized  stopband  frequency  is  obtained,  (in  this  case  0.757r), 
the  design  chart  is  used  to  determine  the  lowest  order  filter  that  will  give: 


3dB  at  9C  =  0.5tt 
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Turning  to  Figure  2.8,  it  can  be  seen  that  the  curve  for  N  =  2  meets 
these  requirements.  At  this  point,  it  is  worth  noting  that  a  second-order  Butter- 
worth  filter  was  also  determined  to  be  required  using  the  traditional  approach  in 
the  previous  section. 

In  Table  2.7  it  can  be  seen  the  prototype  low-pass  filter  transfer  func- 
tion for  N  =  2  is: 

LPpW       3.4l22+0.59       3.4l22+0.59  y"'     } 

Step  5:  Once  the  prototype  filter  is  obtained  from  the  design  curves,  it  is  neces- 
sary to  determine  the  actual  transfer  function  that  will  meet  the  design 
specification,  such  that  0'c  is  0.27T  and  9'a  is  0.47T. 

Using  the  digital  filter  frequency  transformations  of  Table  2.5,  the 
transfer  function  is  determined  as  follows: 

HLP{z)  =  HLPp{z) 

z2  +2z  +  \ 


ZA\z2  +  0.59 

z2  +  2z  +  1 


(2.35) 


UMz2  -17.02  +  6.14 
which  is  the  same  transfer  function  obtained  in  the  previous  section. 
Step  6:     A  computer  generated  frequency  response  plot  of  the  filter  is  used  to 

verify  that  the  design  fulfills  the  design  specifications;  Figure  2.9  reveals 

that  this  is  indeed  the  case. 

D.    LOWPASS  PROTOTYPE  TO  HIGHPASS,  BANDPASS  BAND- 
STOP  FILTERS 

Once  the  lowpass  prototype  filter  is  obtained,  a  highpass,  bandpass  or  band- 
stop  filter  can  be  derived  using  the  digital  frequency  transformations  of  Table  2.5. 
In  the  next  section,  an  example  is  presented  that  illustrates  this  process. 
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1.    Chebyshev  Filters  -  A  Direct  Design  Approach 

As  stated  in  the  introduction  to  this  section,  this  design  approach  is  of 
course  applicable  to  Chebyshev  filter  design  problems.  Again  using  the  bilinear 
transformation  techniques  cited  in  [2],  the  digital  filter  coefficients  corresponding 
to  the  analog  filter  coefficients  for  prototype  Chebyshev  filters  were  determined. 
The  frequency  responses  of  the  resulting  transfer  functions  listed  in  Table  2.8  were 
then  plotted  to  obtain  a  set  of  design  curves  with  a  cutoff  frequency  of  7r/2  (see 
Figures  2.10  thru  2.12). 

The  design  problem  of  the  previous  section  will  be  used  to  illustrate  ap- 
plication of  this  new  technique  (Example  2.2). 

The  problem  statement  called  for  a  Chebyshev  digital  bandpass  filter  to 
be  designed  to  meet  the  following  specification: 

•  1  dB  ripple  in  the  range  of  600  to  900  Hz 

•  Sampling  frequency  is  3  kHz 

•  Maximum  gain  of  -40dB  for  0  <  /  <  200  Hz 

Step  1:    Convert  the  critical  analog  design  frequencies  to  digital: 

Sampling  frequency:  fs  =  3  kHz 

Lower  Passband  frequency:  ft  =  600  Hz 

Upper  Passband  frequency:  /„  =  900  Hz 

Stopband  frequency:  fa  =  200Hz 

d\   =  weT  =  wt/fs  =  2(60°)7r  =  0.4tt  =  1.26  rad 
1  u  3000 

0'u  =  wuT  =  wu/f3  =  2(39Q00°0)7r  =  0.6tt  =  1.88  rad 

9(900W 

S'a  =  waT  =  walfs  =     30Qq     =  0.1337T  =  0.418  rad 
Ripple  band  center  frequency:  6'0  =  y/0'e8'u  =  \/2.35  =  1.54  rad 
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TABLE  2.8 

LOWPASS  PROTOTYPE  CHEBYSHEV  FILTER  TRANSFER  FUNCTIONS 

(due  to  reference  4) 
1/2  -  dB  ripple  (e  =  0.3493) 


N  HW  H{z) 

-i  2.863  2.863z  +  2.863 

3  +  2.863  3.8632  +  1.863 

9  1.43 1.43z2+2.86z+1.43 

^  32  +  l.  425s+l.  516  3.941z2  +  1.032z+1.091 

o  0.716 0.716z^+2.148z^  +  2.148z  +  0.716 

°  33  +  1.25332+l. 5353  +  0. 716  4.5  04z3-0  .5  70  z2+2  .360z-0  .56  6 

a  0.944 0.358z4  +  1.432z3+2.148z2+1.432z  +  0.358 

*  34  +  1.19733  +  1.71732  +  1.025s  +  0.379  5.3  18  z4  -2  .82  8z3+4  .840  z2  -2.1  40  z  +  0  .874 


1  -  dB  ripple  (c  =  0.5088) 


TV  H[s)_  H(z) 

i  1.965  1.965z  +  1.965 

3+1.965  2.965z  +  0.965 

o  0.983  0.983z2  +  1.966z  +  0.983 

L  32  +  l.  0983+1. 103  3.201z2+0.206z  +  1.005 

q  0.491 0.4  91  z3  +  l  .4  73  z2  +  l  .473  z  +  0  .49  1 

33+0.98832  + 1.2383  +  0. 49  1  3.7  17  z3 -1  .2  77  z2 +2  .24  7  z -0  .75  9 

a  0.891 0.2  4  6z'1+0.984z3  +  1.4  76z2  +  0.9  84z  +  0.24  6 

34+0.9  53  33  +  l   45432 +0.743  3  +  0.276  4  .4  26  z4 -3  .3  1  6  ;3 +4  .7  48  z2 -2.4  76  z+ 1  .034 


2  -  dB  ripple  (e  =  0.7648) 


N  H(s)  H(z) 

I  1.308  1.308Z  +  1.308 

3+1.308  2.308z  +  0.308 

9  0.506  0.506z2  +  1.012z  +  0.506 

32+0. 8043+0. 637  2.441z2-0.726z  +  0.833 

q  0.327 0. 3  27  z3 +0. 981  z2 +0.981  z  +  0. 327 

33 +0.73  8. s2  +  l. 0223  +  0. 327  3.087z3-1.7  35z2  +  2.221z-0.95  7 

A  0.164 0.164z4+0.656;3+0.9  84z2  +  0.6  56z  +  0.164 

34+0.7  16  33  +  l  .25  63  2 +0.5  17  3  +  0. 20  6  3.6  95  z4 -3  .57  4  z3 +4  .7  24  z2 -2.778  z  +  1  .22  9 
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TABLE  2.8  (Continued) 

LOWPASS  PROTOTYPE  CHEBYSHEV  FILTER  TRANSFER  FUNCTIONS 

1/2  -  dB  ripple  («  =  0.3493) 

K  £(£l  #[£} 

5 


0.179 

'*+' 

)73« 

•  +  1 

937 

.3+1.310. 

'  +  0.7S3.+0.179 

OOP 

.'+1.159. 

'  +  2.172, 

.'  +  1 

.600.*+l. 

172.-+G  432.+0  096 

O.MS 

•  +] 

.151 

(l+2. 

•113/ 

+1 

.869 

i'  +  1.648r 
0.01? 

>+0.7S6.3+0.282, 

p+U 

045 

6 


+0.1S3.+  0.024 


1-dB  ripple  (c  =  0.5088) 
ff(<)  ff(z) 


r'+0.937. 

•  +  1 

689 

.3  +  0.974. 

0  061 

3  +  0.  SSI.  +  0.123 

i'+0.928»'+1.931 

>*+l 

.202»3+0.939»3+0.307.+0.069 

.923, 

-'  +  2.176. 

>  +  l 

.429. 

p«+1.3M. 
0  015 

3  +  0.549.3+0.214, 

> +0.031 

0.031.-T+O.215.-*+O  64S.-'  +  l  075.--  +  1  075- 3  +  0.645.-:' +  0  21S.  +  0  031 
7. 679, :- IS. 286,'  +  26. 242, '-29. 121, '+25. 126, 3-1 5.812,3 
+6.919.-1.816 

015.-'  +  OI2.-'+0  4:.-_'+0  84.-,+  105,<+0.84,3+0.42.-3  +  0.12,+0.01S 
9. 2M:'  -22.255.  •  +41.993.-'  -53. 661  iJ  +53.51  7.-«  -40.61 l? 
+23. 242, '-9.271, +  2. 196 


2  -  dB  ripple  («  =  0.7648) 


5 

008? 

«'+0.707»<  +  l. 500,3  +  0. 694.3  +  0.459.+  0. 082 

6 

0.041 

»'+0.701»»  +  1.746.«+0. 867.^  +  0. 772.3 +0.210.+0. 051 

7 

0.020 

. '+0.698.'+ 1.994. '+1.039. '  +  1.44. J +0.383. 3  +  0. 167.  +  0. 020 

8 

0.010 

0  041, '+0. 345^  +  0  612.--  +  0  816,3  +  0  SH-"3 +0  245.+  0  04  1 
5.347,'-9. 604,'+ 15. 210,' -IS. 074, 3  +  11.297,3-5.677, +  1.790 

0.020,T+0  I43.-'+0  428.-'  +  0  714, '+0  714.-'+0  428, 3+0  143, +  0.020 
6. 445," -14. 242, '  +  2S. 034, '-29. 203, '+26. 062/3-1 7. 08£, 3 +  7. 765. -2. 16S 

0.010.-'+0  082,T+0  266.-'  +  0.571.-'  +  0  714,,  +  0  571.-'+0  286.-3+0  082.+  0.010 
7. 772, '-20. 397,  '  +39. 593,'  -52. 787, '+54. 371,  '-42. 699, •> 
+  25.30, 3  -10. 46S.  +  2. 615 
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Step  2:  Again,  the  design  chart  is  based  on  a  ripple  edge  frequency  of  7r/2,  necessi- 
tating the  normalization  of  the  design  specification  frequencies.  Looking 
at  the  Table  of  Digital  Filter  Frequency  Transformations  (Table  2.5),  the 
following  frequency  transformations  apply  for  a  bandpass  filter: 

Prototype  filter  from  desired  filter: 

j29'a  _  2ak    j0'a     ,     fc^ 

eJ9*  = fc+i   1_*±J_  (o  36) 


where: 

_cos(9'u/2  +  9'e/2) 

a      cos(6'J2-0'e/2) 


k  =  tan(0c/2)cot 


2        2  J 


9 a  is  the  stopband  frequency  used  to  determine  the  filter  order  N . 

9C  is  the  ripple  edge  frequency  on  the  design  chart,  in  this  case  0C  =  7r/2. 


For  this  problem: 

#a   =? 

9'a  =  0.133tt 

9'u  =  0.6tt 

9'e  =  0.4tt 

9C  =  0.5tt 

_  cos(0.3tt  +  0.2tt)  _  cos(0.5tt) 
01  ~  cos(0.37r-0.27r)  ~  cos(O.Itt) 
k  =  tan(7r/4)cot(0.l7r)  =  3.07 
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Since  a  =  0 


e'"2'-  +  j=f  e^«  +  0.509 

1  +  IrreJ2^  "      1  +  (0.509)c>Mi 
.cjo.836  _  0>509       I.394e-J2.58 


l  +  0.509eJ°-836         1.394e>0-282 


(2.37) 


=$>      9a  =  —2.86  rad  =  2.86  rad  (due  to  symmetry) 

Step  4:      Determine  the  lowest  order  Chebyshev  lowpass  filter  prototype  from  the 
design  chart  that  gives: 


3  db  at  9C  =  0.5tt  rad 


less  than   -  40  db  at  6n  =  2.86  rad 


Turning  to  Figure  2.11,  the  filter  order  is  determined  to  be  N  =  3.  Again, 
this  agrees  with  the  order  obtained  using  the  approach  in  the  previous  section. 
From  Table  2.8  the  lowpass  prototype  transfer  function  for  N  =  3  is: 

0.491*3  +  1.473*2  +  1.473*  +  0.491 
HLP'(Z)  =  3.717*3-  1.277*' +  2.247*  -0.759  (2-38) 

Step  5:  As  with  the  Butterworth  design,  the  Digital  Filter  Frequency  Transfor- 
mation Table  2.5,  is  used  to  convert  the  lowpass  prototype  filter  to  a 
bandpass  filter. 


HBp{z)  =  HLPp(z) 


z--       fc,+  i   :  <i+i_ 


i_  2aA  z+  k-l  z2 
k  +  iz+k+iz 

0.49U3  +  1.473z2  +  1.473-g  +  0.491 

3.717z3  -  1.277*2  +  2.247z  -  0.759 

0.012z6  -  0.036z4  +  0.036;2  +  0.012 

z6  +  2.13624  +  1.76822  +  0.540 


(2.39) 


Again,  this  is  nearly  identical  to  the  transfer  function  obtained  using  the 
'traditional  approach" . 
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Step  6:    Verification  of  design  using  a  computer  generated  frequency  response 
plot.  From  Figure  2.13  it  can  be  seen  the  design  specifications  of: 
$'e  =  0.4tt  =  1.257  rad 
ffu  =  0.6tt  =  1.885  rad 
0'a  =  0.133tt  =  0.418  rad 
are  met. 
2.    Elliptic  Filters  -  A  Direct  Design  Approach 

As  was  shown  in  Example  2.3,  the  design  of  elliptic  filters  is  algebraically 
intensive,  an  extremely  undesirable  feature,  especially  for  high-order  filters.  Elliptic 
filters  characteristically  exhibit  a  sharper  cutoff,  that  is,  a  narrower  transition  width 
for  a  given  filter  order  than  their  Butterworth  or  Chebyshev  counterparts.  For  this 
reason  they  are  very  popular,  making  the  pursuit  of  a  more  efficacious  design 
approach  a  worthwhile  endeavor. 

Just  such  an  approach  was  found  by  taking  the  analog  prototype  filters  of 
Table  2.4  and  normalizing  them  to  have  a  passband  ripple  edge  (critical)  frequency 
of  one.  They  were  then  transformed  to  digital  versions  using  the  bilinear  transfor- 
mation, resulting  in  digital  prototype  filters  with  a  critical  frequency  of  rc/2  (Table 
2.9). 

Thus,  the  amount  of  calculation  involved  in  the  design  of  a  digital  elliptic 
filter  is  cut  in  half  because  the  requirement  for  finding  an  analog  prototype  is 
eliminated.  The  designer  need  only  convert  the  filter  design  specifications  to  digital 
frequencies,  find  the  digital  lowpass  prototype  filter  that  meets  these  requirements 
(from  Table  2.9),  and  then  find  the  actual  filter  transfer  function  through  the  use 
of  the  digital  frequency  transformations  of  Table  2.5. 
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Two  examples  will  be  presented  that  illustrate  this  process;  however,  a 
more  detailed  explanation  of  how  the  transfer  functions  of  Table  2.9  were  arrived 
at  is  in  order. 

For  purposes  of  explanation,  a  second-order  filter  with  0.5  dB  passband 
ripple,  and  a  stopband  gain  of  —20  dB  will  used. 

Turning  to  Table  2.4,  the  lowpass  prototype  transfer  function  is  of  the 

form: 

H0s2  +H0Aqi 


Hlpp(s)  = 


s2  +  Bus  +  Bqi 


where, 


Ho  =0.100220 
A0i  =  5.33789 
Boi  =  0.566660 
Bn  =  0.809390 


therefore, 


Hlpp{s) 


0.1s2 +0.5338 


(2.40) 


s2  +  0.8094.S  +  0.5667 

This  prototype  has  a  value  for  R  of  2.76261.  The  parameter  R  is  the  ratio 
of  the  stopband  frequency,  W2P,  to  the  passband  cutoff  frequency,  u;lp(i.e.,  R  = 
W2P/  wiP)-  Thus,  it  describes  the  sharpness  of  the  transition  region.  A  large  value 
of  R  is  indicative  of  a  broad  transition  region,  while,  conversely,  a  small  value 
corresponds  to  a  narrow  transition  region.  As  expected,  the  higher  the  filter  order, 
the  narrower  the  transition  region,  and  the  smaller  the  value  for  R.  Figure  2.14 
(due  to  [5]),  illustrates  the  magnitude  squared  frequency  response  of  the  normalized 
lowpass  elliptic  filters  of  Table  2.4. 

For  this  particular  example,  w\p  and  wiP  are  subject  to  the  constraint 
that  their  ratio  R  be  2.76261.    Furthermore,  it  should  also  be  noted,  the  filters 
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Figure  2.14.      Magnitude  Squared  Frequency  Response 
of  a  Normalized  Lowpass  Elliptic  Filter 

comprising  this  table  have  been  normalized,  so  that  the  geometric  mean  of  w\p 
and  W2P  is  one. 

(wlpw2p)  =  l  (2.41) 

Combining  these  two  constraints,  the  following  relationships  between  w\p  and  u>2P 
apply: 

wlp  =  l/\/R  =  0.6016  rad/s  (2.42) 

w2p  =  VR  =  1.6620  rad/s  (2.43) 

The  goal  of  the  direct  design  procedure  is  to  find  a  digital  parallel  to  the 
analog  prototype  filters,  and  tabulate  a  table  of  digital  lowpass  prototypes. 

Step  1  of  this  process  is  to  normalize  the  transfer  functions  of  Table  2.4 

so  that  they  all  have  a  value  for  w\p  of  one.  This  is  depicted  Figure  2.15. 

For  this  example: 

0.132  4- 0.5338 
Hlp'<9)  =  ,'+0.8094, +  0.5667  (2'44) 
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Figure  2.15.      Normalization  of  Elliptic  Filter 


where 


wlp  =  l/y/R  =  1/^2.76261  =  0.6061 

To  normalize  Hlpp(s)  to  w\p  =  1,  let  s  =  s/yR  =  0.6061s 

0.1(0.6061s)2  +0.5338 


Hlpp{s) 


(0.6061s)2  +  0.8094(0.60616)  +  0.5667 
0.0367s2  +0.5338 


(2.45; 


0.3674s2  +  0.4906s  +  0.5667 
where  now,  w\p  =  1,  and  w2p  =  R,  since  R  =  W2P/wip. 

The  normalized  analog  prototype  was  then  converted  to  a  digital  prototype 

through  the  use  of  bilinear  transformation. 


Hlpf{z)  =  H[ypp(s) 


0.0367  (ff^-)    +0.5338 


0.367  (7+^)    +  0.4906  (±=±)  +  0.5667 
1.575z2  -2.752  +  1.575 


(2.46) 


3.9U2  +  1.132  +  1.22 
The  digital  passband  critical  frequency  is  8lp    =   7r/2,  since  an  analog 

frequency  of  w\p    —   1   corresponds  to  digital  frequency  of  7r/2  when  using  the 

bilinear  transformation. 
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Table  2.9  contains  the  compiled  results  of  applying  this  process  to  some  of 
the  analog  prototype  transfer  functions  of  Table  2.4.  As  with  analog  elliptic  filters, 
special  relationships  exist  between  the  passband  critical  frequency,  9\p,  and  the 
stopband  frequency,  02P,  of  digital  elliptic  filters. 

Using  the  digital  frequency  transformations  of  Table  2.5,  the  following 
values  for  B\p  and  $2P  are  found  based  on  the  analog  prototype  values  for  W\p  and 
w2p- 


For9lp, 


For  0 


2P, 


-jl  +  1 

n/4 

--  le^/2 

*  Blp  =  tt/2 

i02p 

jR  +  l 
-jR  +  1       le 

Lejtan- 
-j tan-1 

1  R 

[(-R) 

=  le>2t 

02P 

=  2  tan-1  R 

(2.47) 


(2.48) 


Again,  R  is  the  transition  region  parameter  described  by  the  ratio  of  w2p 
to  Wip  of  the  analog  version  of  the  filter.    Table  2.9  presents  a  summary  of  the 
values  for  9\p,  $2P,  and  R  for  some  of  the  digital  lowpass  prototype  elliptic  filters. 
Example  2.6 

Once  again  the  design  example  of  the  previous  section  will  be  used  to 
illustrate  the  new  direct  design  procedure.  The  problem  statement  is: 

Design  an  elliptic  lowpass  filter  that  meets  the  following  specifications: 

•  passband  ripple  of  0.5  dB 

•  passband  ripple-edge  frequency  of  2  kHz 

•  stopband  gain  of  at  most  —20  dB  for  frequencies  greater  than  6  kHz 

•  sampling  frequency  is  20  kHz. 


76 


TABLE  2.9 
ELLIPTIC  LOWPASS  PROTOTYPE  FILTERS 

Passband  ripple  =  0.5  dB  ;  Stopband  gain  =  —20  dB 
(due  to  reference  5) 


N       9ip_         dj^  R 

2  1.5708  2.447  2.7626 

3  1.5708  1.9157  1.9157 

4  1.5708  1.6145  1.0447 

5  1.5708  1.5862  1.0155 


Transfer  Functions 


N  H(z) 

o  1.575z2  +  2.75s  +  1.575 

^  3.91z2  +  1.13z  +  1.22 

o  1.276z3  +  2.367z2+2.376z  +  1.276 

°  4.626z3+0.008z2+3.030z-0.352 

a  1.4179z4+2.366z3+3.4362z2+2.336z  +  1.4179 

*  5.885z4-1.0958z3  +  6.5858z2-1.0954z  +  1.338 

r  1.794z5+2.9512z4+4.8532z3+4.8532z2+2.9512z  +  1.794 

°  7.9516z5-2.314z4+12.7192z3-3.1856z2+4.9276z-0.902 
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Step  1:    Find  the  critical  digital  design  frequencies. 

=  2^2,(2x103) 
lp         fa  20  x  103 

2^      2,(6  x  1Q3)  = 
2p  /,  20  x  103 

Step  2:    To  determine  the  filter  order  of  the  prototype  filter,  find  a  and  92p ,  which 

are  defined  as  follows  from  Table  2.5. 


sin  (flc/2- 0^/2) 
sin(0c/2  +  0'    /2) 


where  0C  is  the  critical  frequency  of  the  prototype  filter  i.e.,  9C  =  tt/2. 

Therefore, 

sin(0.57r/2-0.27r/2) 
a 


sin(0.57r/2  +  0.27r/2) 
_  sin  (0.7854 -0.31416)  _  0.45399  (2.49) 

~  sin  (0.7854  +  0.31416)  ~  1.09956 
=  0.5095 
Find  92p ,  the  stopband  frequency  for  the  lowpass  prototype  filter. 


e>9*p  =   € 


1  -  aeje*p 
eji-855  _  0,5095 

"  1  -  (0.5095)e>1-855  (2.50) 

_       i-^QQe _  -.J2.6785 

l^SSe-^0-3965 
=>92p  =  2.6785  rad 

Step  3:  From  the  design  chart  (Table  2.9),  determine  the  lowest  order  elliptic 
lowpass  filter  prototype.  Table  2.9,  indicates  that  a  second-order  fil- 
ter has  a  stopband  frequency  92p  of  2.447  rad,  which  meets  the  design 
specifications.  Thus,  the  lowpass  prototype  transfer  function  is: 

1.575;2+ 2.75; +  1.575 
HlPp{z)  ~    3.91*'  + 1.13* +  1.22  (2-51) 
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Step  4:    Find  the  actual  filter  transfer  function,  using  the  prototype  transfer  func- 
tion and  the  digital  filter  frequency  transformations  of  Table  2.5. 


Hlp(z)  =  Hlpp{z 

1.57522  +  2.752  +  1.575 


3.9U2  +1.13* +  1.22 


(2.52) 


_  0.5829z2  +  0.2541*  +  0.5829 
3.65U2  -  3.80422  +  1.6593 

or 

TT     ,  N       0.1596z2  +  0.0696*  +  0.1596  ,n  roN 

HLP{Z)  =         z*-  1.042.  +  0.4545  (2'53) 

The  transfer  function  of  Equation  (2.52)   is  identical  to  the  transfer 

function  obtained  using  the  traditional  design  technique  (see  previous 

section). 
Step  5:    Verify  the  design  by  obtaining  a  computer  generated  frequency  response 

plot  (see  Figure  2.16). 
To  further  illustrate  this  direct  design  method,   an  additional  example 
follows  involving  the  design  of  a  bandpass  filter. 
Example  2.7 

A  digital  elliptic  bandpass  filter  is  to  be  designed  to  meet  the  following 
specifications: 

•  0.5  dB  ripple  in  the  frequency  range  600  Hz  <  /  <  900  Hz, 

•  sampling  frequency  of  3  kHz,  and 

•  maximum  gain  of  -20  dB  for  0  <  /  <  200  Hz. 
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Figure  2.16.     Frequency  Response  for  Elliptic  Lowpass  Filter 


SO 


Step  1:      Find  the  critical  digital  frequencies  from  the  design  specifications: 
1        f3  3000 

,a  =  ^  =  M9oo)         =lg8rad 

fs  3000 

2p         /a  3000 


^  =  /^  =  v/(1.25)  (1.88)  =  1.53  rad 

Step  2:  Again,  since  the  design  chart  is  based  on  a  cutoff  frequency  of  0C  =  ir/2, 
the  design  specifications  need  to  be  normalized  to  determine  the  lowpass 
prototype  filter  from  Table  2.9. 

From  the  Table  of  Digital  Frequency  Transformations  (Table  2.5), 
the  following  frequency  transformations  apply: 

Prototype  filter  from  desired  filter: 

J20'       _   2ak  J9'2  fc-i 

e^r  =  -- y6     -    +  y  (2.53) 

1       fc+ie  +  Jfc+ie 

where: 

_cos(^/2  +  ^/2) 

a       cos(0'u/2-^/2) 
fc  =  tail  (0c/2)  cot  (0^2 -0J/2) 
6'2p  is  the  stopband  frequency  used  to  determine  the  filter  order,  JV. 
0C  is  the  cutoff  frequency  of  the  prototype  filter,  8C  =  n/2. 
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For  this  problem: 

6'2p  =  0.133tt 

0'u   =0.6tt 

9'e   =0.4tt 

6C  =  0.5tt 

_  cos(0.37r  +  0.27r)  _  cos(0.5)7r 
a  ~  cos(0.3tt  -  0.2tt)  ~  cos(O.Itt) 
a  =  0 

fc  =  tan(0.57r/2)cot(0.67r/2  -  0.4tt/2)  =  3.07 

Step  3:    Find  92p  to  determine  the  lowest  order  prototype  filter  that  may  be  used. 
Since  a  =  0 

JH9-      ^^  +  fe*  eJ^  +0.509 


1  +  f^e'2'2*  1  +  (0.509)ei2^p 

eJo.836+0509  i.394e->258       1    _,286         (2'55) 

—  =le   J  " 


1  +  (0.509)eJ°-836        1.394e>0-282 
==►  02p  =  2.86  rad 

Step  4:    From  Table  2.9,  determine  the  filter  order  of  the  prototype  filter  that 
corresponds  to  the  above  value  of  #2P-  It  can  be  seen  from  the  table  that 
a  second-order  filter  fulfills  the  design  specifications. 
Again,  the  lowp?,ss  prototype  transfer  function  is: 

1.575z2+2.75z  + 1.575 
*"*W  =    3.91,'  + 1.13, +  1.22  (2-56) 
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Step  5:  As  in  previous  design  problems,  the  prototype  filter  transfer  function  is 
converted  to  the  actual  filter  transfer  function  through  the  use  of  the 
digital  frequency  transformations  of  Table  2.5. 


HBP(z)  =  HLPp(z 


z=-- fc  +  i   r  Eu 

1.575*2  +2.752  +  1.575 


3.9122  +  1.132  +  1.22 
0.5833z4  -  0.2558-z2  +  0.5833 


(2.57) 


-a     ?2+Q5Q9 

-0.509z2-l 


3.650924  +  3.7996z2  +  1.6579 
Step  6:     Obtain  the  frequency  response  plot  for  design  verification  (see  Figure 
2.17). 

In  summary  then,  this  chapter  has  presented  three  commonly  used  recursive 
filter  types:  Butterworth,  Chebyshev  and  elliptic.  The  type  chosen  by  the  designer 
is  dependent  on  the  requirements  of  his/her  particular  design  problem. 

Butterworth  filters  exhibit  a  flat  frequency  response  characteristic,  but  do  not 
have  as  steep  a  transition  region  for  the  same  order  filter  as  do  Chebyshev  and 
elliptic  filters.  With  Chebyshev  and  elliptic  filters  the  steepness  of  the  transition 
region  is  attained  by  accepting  a  degree  of  ripple  in  the  passband. 

All  three  filter  types  may  be  designed  using  either  the  traditional  approach, 
involving  an  analog  prototype  filter  transfer  function,  or  the  direct  design  approach 
that  uses  digital  prototype  filter  transfer  functions. 
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Figure  2.17.     Frequency  Response  for 
Elliptic  Bandpass  Filter 
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III.    NONRECURSIVE  FILTER  DESIGN 

A.    INTRODUCTION 

The  nonrecursive  realization  of  digital  filters  is  desirable  because  it  can  result 
in  two  very  attractive  features:  linear  phase  and  the  absence  of  stability  problems 
because  all  of  the  poles  are  at  2  =  0. 

Nonrecursive  filters  are  normally  characterized  by  a  finite  duration  impulse 
response;  consequently,  historical  methods  of  design  involve  the  analytical  deter- 
mination of  the  filter  coefficients  by  expanding  the  desired  frequency  response  in  a 
Fourier  series  and  then  truncating  the  series  to  the  desired  filter  length  (order).  A 
disadvantage  inherent  in  this  design  method  is  the  presence  of  an  overshoot  that 
occurs  at  discontinuities  in  the  desired  frequency  response  (known  as  Gibbs'  phe- 
nomenon). The  use  of  window  functions  was  devised  as  a  remedy  to  this  problem; 
examples  are  the  Kaiser,  Hamming  and  von  Hann  windows  [14]. 

Another  popular  method  of  FIR  filter  design  is  frequency  sampling,  whereby 
the  desired  filter  frequency  response  is  sampled,  and  then  the  inverse  discrete 
Fourier  transform  of  these  sample  values  is  determined  to  find  the  filter  impulse 
response. 

In  this  chapter  both  of  these  methods  will  be  presented,  and  several  design 
examples  given,  but  first  a  review  of  the  theoretical  background  of  nonrecursive 
digital  filter  design  by  Fourier  methods  is  in  order. 
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B.    BACKGROUND 

The  design  of  nonrecursive  digital  filters  hinges  on  the  following  relationships: 
A  causal  nonrecursive  system  can  be  described  by  a  difference  equation  of  the 
form  [1] 

k=L 

y(n)=  Y,hkx(n-k).  (3.1) 

jfc=0 

For  an  input  x(n)  =  ejn*  the  steady  -  state  system  output  is 

y99(n)  =  jn9H{J9)  (3.2) 

where  H(eJ  )  is  the  system's  frequency  response. 

Expanding  the  right  side  of  Equation  (3.1)  for  the  same  input  x(n)  =  e^n6 , 
yields  the  following  steady  -  state  output 

y„(n)  =  £  fee*-**  =  £  h^'e-^ 


(3.3) 


Jt=0  Jfc=0 

=  b0eine  +  bxe?net->9  +  b2e>n$e-32e  +  ....  +  bLe>n9e-'L9 

=  e?n6  [b0  +  he''9  +  fee"'""  +  . . .  +  bLe''L$]  . 
Equating  Equations  (3.2)  and  (3.3)  gives 

e>n9H{e>9)  =  e>n9  [b0  +  6lC"^  +  b2e~]29  +  ...  +  bLe^L0]  (3.4) 

which  implies 

H(e>9)  =  Y,  Ke-Jn6.  (3.5) 

n=0 

But,  by  definition  the  frequency  response  of  an  LTI  system  is 

n=oo 

H(e*9)  =   £  A(n)c->»*.  (3.6) 

n=oo 

For  a  causal  filter  (h(n)  =  0,n  <  0),  with  a  finite  number  of  delays  (n  =  L), 
however,  the  frequency  response  definition  is 

n=L 

H{e>9)  =  £  h(n)e-J"9.  (3.7) 
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Comparing  Equations  (3.5)  and  (3.7)  produces  the  following  relationship  that  is 
the  crux  of  nonrecursive  filter  design 

h(n)  =  bn.  (3.8) 

Equation  (3.8)  establishes  a  direct  relationship  between  the  unit  impulse  response 
of  an  FIR  filter  and  the  coefficients,  6n,  in  the  system  difference  equation. 

The  goal  of  the  design  techniques  presented  in  this  chapter  is  to  determine  the 
filter  coefficients  (or  weights),  &o,  &i , . . . ,  &£,  for  a  desired  frequency  response  charac- 
teristic, Ho(^9)-  These  techniques  include:  Fourier  Coefficient  Design  for  lowpass, 
highpass,  bandpass  and  bandstop  filters;  Windowing;  and  Frequency  Sampling. 

C.    FOURIER  COEFFICIENT  DESIGN 

As  stated  in  the  introduction  we  need  to  determine  the  filter  coefficients,  which, 
in  the  case  of  nonrecursive  filters,  are  also  the  coefficients  of  the  unit  impulse 
response,  h(n).  Thus,  a  relationship  between  the  desired  frequency  response  and 
the  impulse  response  needs  to  be  established. 

Beginning  with  the  desired  frequency  response 

n=oo 

HD(eje)=    Y,    h(")e~jne-  (3-9) 

n=— oo 

It  can  be  shown,  [1],  that  h(n)  of  Equation  (3.9)  may  be  written, 

JOo 


2n+60 

h(n)  =  (1/2tt)  /  HD(ejeyn9d9.  (3.10) 


Equation  (3.9)  is  recognized  as  the  Fourier  series  expansion  of  the  function  Hd(^&)i 
where  the  h(n)  are  the  Fourier  coefficients  (hence  the  source  of  the  name  for  this 
design  technique). 

Depending  on  the  type  of  filter  that  is  desired,  this  expression  can  be  reduced 
to  the  following  forms: 
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(3.n; 


Lowpass  Filters: 

^Lp(rc)  =  (K/irn)sm(n0c) 

K  —   desired  magnitude  in  the  passband 
n  =  0,±l,±2,... 
9C  =    the  desired  cutoff  frequency 

The  number  of  coefficients  is  truncated  to  correspond  to  the  desired  filter  order; 
as  expected,  the  higher  the  order,  the  better  the  approximation  to  the  desired 
frequency  response. 
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0 

Figure  3.1.      Ideal  Lowpass  Filter  Frequency  Response 


Highpass  Filters: 


hHp(n)=hLP(n)(-l)n 


(3.12) 
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Figure  3.2.     Ideal  Highpass  Filter  Frequency  Response 

Bandpass  Filters: 

hBp(n)  =  [2cos(n0o)]  hLP(n)  (3.13) 

where, 


Bu-Si 

2 

6u+0t 


>o  = 


;      #u  =    upper  passband  frequency 
;      0(  =   lower  passband  frequency 


k    -e.-e. -a, 


e-j    °0     °-      K      y 


Figure  3.3.      Ideal  Bandpass  Filter  Frequency  Response 

Bandstop  Filters: 

fcz?s(0)  =  K  -  hBP(0) 

hBS(n)  =  -hBP(n);n  =  ±1,±2,... 


(3.14) 
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HD(e'     ) 


Figure  3.4.     Ideal  Bandstop  Filter  Frequency  Response 

As  can  be  seen  in  the  previous  impulse  response  expressions  for  the  various 
filter  types,  they  are  all  based  on  a  lowpass  filter  prototype,  hLP(n).  Thus,  the 
design  steps  for  the  ideal  lowpass  prototype  filter  can  be  summarized  as  follows, 
bearing  in  mind  that  once  the  lowpass  filter  coefficients  are  determined  they  may 
be  transformed  to  the  coefficients  for  a  highpass,  bandpass,  or  bandstop  filter,  if 
desired. 

1.    Lowpass  Prototype  Filter  Design  Procedure 
SteP  1:    Translate  the  design  specifications  to  those  of  a  lowpass  prototype.  De- 
termine the  desired  critical  frequency,  6C,  and  the  passband  magnitude, 
A'. 

Step_2:     Find  the  lowpass  filter  coefficients  given  by  Equation  (3.S),  that  is, 

hLP(n)  =  bn  =(K/Tm)sin(nOc)i    n  =  0,  ±1,  ±2,. . . ,  ±L         (3.15) 

with  L  =  (N  -  l)/2. 
Step_3:    Shift  hLP(n)  to  the  right  by  L  terms  to  mate  the  filter  causal. 


hLP{n)  =  [A'/7r(n  -  L)\  sin  [(a  -  L)9C  J ;     n  =  0, 1, 2, . . . , 2L       (3. 


16) 
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Step  4:  Transform  the  coefficients,  hLp(n),  to  lowpass,  highpass,  bandpass  or 
bandstop,  as  desired.  Implement  the  design,  and  compare  the  frequency 
response  obtained  with  the  original  specifications  to  ensure  that  these 
specifications  are  met. 

As  an  example,  suppose  a  highpass  filter  with  a  passband  of  unity  gain  for 
frequencies  greater  than  10  kHz  is  desired.  The  system  sampling  frequency  is  50 
kHz. 

Step  1:     Cutoff  frequency,  9'c  =  wT  =  ^  =  fffipj  =  0.4tt 

Passband  gain,  K  =  1.0 

Translating  to  lowpass  prototype  specifications: 

&c  =  7r  —  9C  ;     where  9C  is  the  cutoff  frequency  of  the 
lowpass  prototype  and  9'    is  the  cutoff 
frequency  of  the  highpass  filter 
=>  9r         =  7T  -  9'    =  7t  -  0.4tt  =  0.6k 

"  c 

Step  2:     From  Equation  (3.11),  the  lowpass  prototype  filter  coefficients  are: 

hLP(n)  =  sin(0.67m)/(7rn);  n  =  0,  ±1,  ±2, . . . ,  ±10 

for  a  21 -coefficient  filter,  while  the  highpass  filter  coefficients  are  obtained 
from  Equation  (3.12).  The  resulting  frequency  response  is  illustrated  in 
Figure  3.5. 
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D.    WINDOWS 

It  has  been  shown  that  the  expression  for  the  desired  frequency  response, 
HD{z*e),  can  be  written  in  terms  of  a  Fourier  series.  Equations  (3.9)  and  (3.10) 
are  rewritten  below  for  convenience. 

n=oo 

HD(ej0)=    J2    h(n)e-jn6d0  (3.9) 

n=—  oo 

where 

27T+0O 

h(n)  =  (1/2tt)  /  HD(eje)ejned8.  (3.10) 

JOo 

In  actual  implementation,  however,  the  infinite  sum  in  Equation  (3.9)  must  neces- 
sarily be  truncated,  since  a  filter  cannot  have  an  infinite  number  of  delays.  Also, 
at  points  of  discontinuity  in  the  ideal  frequency  response  (where  the  magnitude 
abruptly  changes  from  1.0  to  0.0,  or  vice  versa),  the  Fourier  series  approximation 
cannot  match  the  desired  frequency  response  exactly,  even  if  an  infinite  number  of 
terms  were  possible.     There  is  always  an  overshoot  of  about  nine  percent  near 
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Figure  3.5.     Frequency  Response  for  Fourier  Coefficient  Highpass  Filter 
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discontinuities.  This  overshoot  is  the  well  known  Gibbs'  phenomenon,  and  is  alle- 
viated through  the  use  of  window  functions.  The  filter  weights,  h(n),  are  modified 
using  one  of  several  available  analytical  expressions  (Hamming,  Von  Hann,  etc.) 
[14]. 

h(n)  =  h(n)  •  w(n)  (3.17) 

h(n)  =  modified  unit  sample  response 

h(n)  =  original  unit  sample  response 

w(n)  —  window  function 

The  modified  value  of  the  desired  frequency  response  is  therefore 

n=L 


n=-L 
n=L 

=    J2    h(n)  •  w(n)e-jn 


(3.18) 


n=-L 


again  where, 


/  iWV 


h(n)  =  (1/2tt)   /     HD(ej6)ejndd8 


which  exhibits  a  gradual  roll-off,  rather  than  the  steep  slope  characteristic  of  the 
ideal  frequency  response.  Thus,  the  tradeoff  involved  when  using  windows  in 
Fourier  design  is  that  a  reduction  in  the  overshoot  caused  by  Gibbs'  phenomenon, 
is  achieved  at  the  expense  of  decreasing  the  sharpness  of  the  frequency  response 
cutoff.  That  is  to  say,  the  passband  and  stopband  ripples  are  suppressed  at  the 
expense  of  a  wider  transition  from  passband  to  stopband;  specifically,  there  is  a 
less  steep  transition.  Figures  3.6  and  3.7  illustrate  these  points;  they  depict  an  un- 
windowed  and  windowed  lowpass  filter  design,  respectively.  For  convenience,  some 
of  the  more  popular  window  functions  are  summarized  [14]: 
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Figure  3.6.       Unwindowed  Lowpass  Filter  Frequency  Response 
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Figure  3.7.     Windowed  Lowpass  Filter  Frequency  Response 
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•  Hamming  Window 

w(n)  =  0.54  +  0.46cos(n7r/I) 
n  =  6,±l,±2,...,±Z 

•  Von  Hann  Window 

w(n)  =  0.50  +  0.50cos(n7r/L) 
n  =  0,±l,±2,...,±L 


(3.19) 


(3.20) 


•  Kaiser  Window 


mi-(2n/L>W» 


hP  (3.21) 

n  =  0,±l,±2,...,±L 

where  Iq{-)  is  the  modified  zeroth-order  Bessel  function,  and  j3  is  a  constant  that 
specifies  a  frequency  response  tradeoff  between  the  peak  height  of  the  sidelobe 
ripples,  and  the  width  of  the  main  lobe. 

E.    DESIGN  OF  A  DIFFERENTIATOR 

The  Fourier  series  design  procedure  is  easily  applied  to  the  design  of  a  differ- 
entiator, which  is  often  used  in  signal  processing  applications  to  track  the  rate  of 
change  of  a  signal.  The  desired  analog  frequency  response  for  a  differentiator  is 

HD(jw)  =  jw.  (3.22) 

Its  counterpart  in  the  digital  frequency  domain  is 

HD(eje)  =  j9/T;9  =  wT.  (3.23) 

Figure  3.8  illustrates  the  ideal  magnitude  and  phase  characteristic. 
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Figure  3.8.     Ideal  Differentiator  Frequency  Response  Characteristic 

The  Fourier  series  design  approach  requires  the  determination  of  the  Fourier 
series  that  approximates  this  ideal  frequency  response  characteristic. 

Substituting  the  desired  frequency  response  (Eqn.  3.23)  into  Equation  (3.10) 
yields 

h{n)  =  (l/2ir)  /    (je/T)ejn6d6 

J-\  (3.24) 

=  (J/2ttT)  I    6ejn6d9. 

Evaluating  this  integral,  the  following  expression  is  obtained  for  the  unit  sample 
response  of  a  differentiator. 

fe(n)  =  (-l)V(nD,n  =  ±l,±2,... 

h(n)  =  0         ,n  =  0 


(3.25) 
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Figure  3.9  depicts  the  frequency  response  of  a  20th-order  differentiator  whose  co- 
efficients for  T  =  1  are  as  follows: 

n        h(n)  n        h(n) 
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Also,  included  is  a  20th-order  differentiator  with  a  Hamming  window,  which,  as 
can  be  seen  in  Figure  3.10,  is  a  very  good  approximation. 

In  summary  then,  the  Fourier  series  design  technique  is  very  effective,  and 
can  be  used  in  many  applications  such  as  the  design  of  a  differentiator  illustrated 
here.  The  use  of  window  functions  is  necessary,  however,  to  reduce  the  Gibbs' 
phenomenon  effect,  which  introduces  a  sacrifice  in  the  designed  filter's  transition 
region. 

F.    FREQUENCY  SAMPLING 

The  window  design  technique  introduced  in  the  previous  section  has  a  distinct 
disadvantage  in  that  the  computation  of  Fourier  series  coefficients  for  the  desired 
frequency  response  can  be  impossible  when  an  analytical  expression  for  the  desired 
filter's  frequency  response  is  unknown  or  extremely  difficult  to  determine. 

Frequency  sampling  is  a  technique  originally  proposed  by  Gold  and  Jordan  and 
later  developed  by  Rabiner  et  al.  as  a  solution  to  this  problem  [9].  This  method 
has  two  distinct  advantages: 
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Figure  3.9.     Unwindowed  21  Coefficient  Bandpass  Differentiator 
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Figure  3.10.     Windowed  21  Coefficient  Bandpass  Differentiator 
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•  It  is  possible  to  design  a  filter  that  approximates  any  desired  frequency  domain 
specifications,  without  ever  having  to  determine  an  analytical  expression  for 
the  filter  frequency  response,  H(e]9). 

•  The  filter  coefficients,  h(n),  can  be  determined  using  the  IDFT  algorithm. 
Before  proceeding  with  an  explicit  summary  of  the  steps  involved,  and  a  de- 
tailed example,  a  discussion  of  the  theoretical  background  will  be  presented. 

In  many  filter  applications  a  sharp  cutoff  amplitude  characteristic  and  linear 
phase  are  desired.  To  ensure  that  these  objectives  are  met  when  designing  a  filter, 
requires  consideration  of  the  number  N  and  location  of  the  equispaced  frequency 
samples.  Discussion  of  these  parameters  relative  to  their  impact  on  the  filter's 
frequency  response  involves  consideration  of  four  separate  cases  [9]. 
Case  1  N  odd,  frequency  samples  at 

0,  =27^/^  =  0,1,2,...,^-!. 


N   =»    9 


Figure  3.11a.      Frequency  Samples  for  N  Odd 
Case  2  N  even,  frequency  samples  at 

0k  =2rr^/iV,;-  =  0,l,2,...,.V-l. 
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Figure  3.11b.      Frequency  Samples  for  N  Even 

Case  3  N  odd,  frequency  samples  at 


6t  .^ii+^i.t.  0,1,2,.  ..,n-i. 


N 


N    =    9 


Figure  3.11c.      Frequency  Samples  for  N  Odd 

Case  4  N  even,  frequency  samples  at 


(A; +  1/2) 

=:    fir 


N 


fc  =  0,l,2,...,iV-l, 
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<f  2^-8/2 


N      = 


Figure  3. lid.      Frequency  Samples  for  N  Even 

It  should  be  noted  that  Cases  3  and  4  do  not  readily  lend  themselves  to  inverse 
DFT  computation  using  the  FFT  algorithm,  since  the  first  frequency  sample  is  not 
at  0  Hz.  For  this  reason  a  detailed  discussion  of  these  two  cases  will  be  omitted, 
however,  information  regarding  these  cases  can  be  found  in  [9].  A  discussion  of 
frequency  sampling  considerations  for  Cases  1  and  2  follows. 
Cases  1  &  2:  Sampling  Theorem 

Recall  that  for  an  FIR  filter  with  impulse  response,  h(n)  —  {/*(0),  h(l), . . . , 
h(N  —  1)}  ,  the  transfer  function  is 

JV-l 


H(=)  =  y,  h^= 


(3.26) 


As  discussed  earlier,  the  impulse  response,  h(n),  can  also  be  represented  in  terms 
of  its  Discrete  Fourier  Transform  (DFT),  since  it  is  of  finite  duration  [9] 

N-l 


h(n)  =  (1/JV)  Y  Hkei2*kn 


/N 


(3.27) 


where     Hk  —  H{z 


;=ei'«*/" 


104 


Substituting  the  expression  for  h(n)  into  Equation  (3.26)  gives  [9] 


N-l 


*<«)-£ 


n=Q    I 


N-l 


(1/N)  J2  Hkej2nkn 


/N 


N-l 


(i/ff)  E 


Ar=0 

N-l 


2jj*. 


,j2irkn/N 


jk=0 
N-l  N-l 


=  (1/N)  J2  E  HkeJ2irkn'Nz-n 

jt=0    n=0 
N-l  N-l 

=  (1/iV)  Y,  Hk  Y,  ej2"kn'Nz-n 


k=0  n=0 

Applying  the  finite  geometric  sum  property  to  the  inner  summation 

N 


N-l 


V^  eJ2nkn/N z-n  _ 


1- 


J2rk/N 


n=0 

N-l 


1- 


e;2xfc/N 


E eJ'2' 


kn/N-n  _ 


1-     * 


Tl=0 


ej2xk/N 


1  _  eJ2*kz-N 


1  _  ej2nk/Nz-l 

Substituting  back  into  Equation  (3.28)  yields  [9] 

N-l 


;ei2,r*  =  l,ifc  =  0,l,2,... 


jfc=0 


(I-*-*)" 


iV 


Jt=o 


ej2nk/N2- 


Let  2  =  e-7^,  to  determine  the  interpolated  frequency  response 

,-jN0/2  („jN0/2  _  jj 


tf(e*)  = 


,-(eJN9/2_e-jne/2j 


N 


jn6/2\   N-l 

f~L  1 


jfc=0 


(ei2*klN)(e-i9) 


_  c-iN»/2   ^J  #fc  (ejN*/2  _  e-J»N/2) 

~~ JV        A*     i-(ci2**/N)(c->«) 


Jfc=0 


(3.28) 
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e-jN6/2   N-l  ffk  ^jN9/2  _  e~jN8/2)  .  e~jwk/N  .  eJ9/2 

=  M  ^  e-jnk/N  .  eJB/2  _  ejnk/N  .  e~j0/2 

Jb=0 

e-jN9/2  .  ejB/2  *zi  Hk  (e>N9'2  -  e->Nel2)  ■  e->vk'N 

=  tf  2^         ej(6/2-*k/N)  _  eJ(0/2-irk/N) 

fc=0 

*  &  *»«-#)'  (3-9) 

Thus,  Equation  (3.29)  expresses  the  interpolated  frequency  response  in  terms  of 
the  sample  values,  Hk,  of  the  desired  frequency  response  [9]. 

In  Case  1,  where  the  number  of  frequency  samples,  TV,  is  odd,  choosing  the 
frequency  samples,  Hk,  to  be  real  and  symmetric  yields  a  real  and  symmetric 
impulse  response 

(N-l)/2 


h( 


»>-f  +  £  ^  «-(£■*)•  (3-30) 

Jk=l  V  7 


The  interpolated  frequency  response,  H  (e-7*),  is  also  purely  real  which,  as  stated 
in  the  beginning  of  this  section,  is  highly  desirable  for  most  filter  applications 

(JV-l)/2 

H  (e>9)  =        Yl       2h  (n) cos  (nd)  ■  (3-31) 

Case  2,  however,  presents  a  problem.  Looking  at  Equation  (3.29),  it  can 
be  seen  that  for  N  even,  the  term  e~i*klN  inside  the  summation  introduces  an 
imaginary  component  into  the  interpolated  frequency  response,  H  (eJ<9). 

For  example,  suppose  the  following  lowpass  filter  frequency  response  is  sampled 
from  —  7r  <  0  <  7r  with  an  even  number  of  frequency  samples,  N  =  4.  From 
Equation  (3.27),  the  following  noncausal  impulse  response  is  obtained 


h(n)  =  (1/4)   V  Hkej2™k/* 

*=-2  (3.32) 


=  (1/4)  [O  +  (l)e-'2^  +  1  +  (l)e1"1** 
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-tc      -0.75k 


HD(e      ) 


0.757C 


Figure  3.12.     Desired  Lowpass  Filter  Frequency  Response 

therefore, 

h{-2)  =  (l/4)[0  +  ey>  +  1  +  e-J>]  =  -0.25 

fc(-l)  =  (l/4)[0  +  ^2  +  l  +  e"^2]  =  0.25 

h(0)  =  (l/4)[0  +  1  +  1  +  1]  =  0.75 

h(l)  =  (l/4)[0  +  e~^2  +  1  +  e^/2]  =  0.25 
and  the  frequency  response  is 

F(e")  =    J]  >>(n)e-» 

n=-2 

=  (-0.25)e-?'2*  +  (0.25)6*  +  (0.75)  +  (0.25)e~* 

=  (-O.25)[cos(20)  +j  sin(20)]  +  O.25[cos(0)  +  jsin(0)]  +  (0.75)      (3'33^ 

+  (O.25)[cos(0)-jsin(0)] 

=  (0.75)  +  (0.5)  cos  6  +  (-O.25)[cos(20)  +  j  sin(20)]. 
The  imaginary  component  introduced  is 

-0.25;  sin(20). 

According  to  Reference  9,  the  amplitude  of  this  component  should  be, 

A  =  (l/N)  £iTt(-l)* 

k=-2 

=  (l/4)[0  -  1  +  1  -  1]  =  -0.25 
which  is  indeed  the  case. 


If  an  odd  number  of  samples,  N  =  5,  were  taken  of  the  same  frequency  re- 
sponse, the  impulse  response  would  be 

2 


h(n)  =  (1/5)  £  Hke**nk'5 


k=-2 

(l/5)[0  +  (l)e-j2^  +  1  +  (l)e'2^  +  0] 


=  (l/5)[l  +  2cos(27rn/5)] 


(3.34) 


(3.35) 


therefore, 

h(0)  =  0.6 

h(±l)  =  (1/5)[1  +  2  cos  2tt/5]  =  0.324 

h(±2)  =  (1/5)[1  +  2cos47r/5]  =  0.124. 
The  frequency  response  is  thus 

H(e*9)  =    £   *(n)e-'»' 

n=-2 

=  (0.324)e'2'  +  (0.124)eJ*  +  0.6  +  (0.124)e"''  +  (0.324)e-'w 
=  (O.324)[cos(20)  +  jsin(20)]  +  (0.124) [cos  0 +  j  sin  0]+ 
O.6  +  (O.124)[cos0- jsinfl]  +  (O.324)[cos(20)  -jsin(20)] 
=  0.6  +  (0.648)  cos  9  +  (0.248)  cos(20) 

and  there  is  no  imaginary  component. 

To  remedy  the  imaginary  component  problem  for  the  even  case,  the  following 
substitution  is  made  for  the  frequency  sample  values  Hk  [9] 

Hk  =  Gkejnk/N.  (3.36) 

The  set  of  Gfc's  is  selected  so  that   G(0)    =    H(0),    G(N/2)    =    0, 
and  G(k)  =  -G(N  -  jb),  k  =  1,  2,  . . .,  {N/2)  -  1. 

The  resultant  impulse  response  is  both  real  and  symmetric,  with  h(n)  =  h 
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(JV  -  1  -  n);  n  =  0,  1,  2,  . . .,  (iV/2)  -  1.  The  values  of  /*(n)  are  given  by, 

w=£+TfH(£)fc+(IH      (337) 

and  the  interpolated  frequency  response  H(e^e)  is, 

(7V-l)/2 

H(ej°)=       Yl       2h(n)cos(n0).  (3.38) 

The  steps  involved,   when  designing  a  filter  using  the  frequency  sampling 
method,  are  as  follows: 

Step  1:  Given  the  continuous  frequency  response  specifications  for  a  desired  filter, 
N  samples  are  taken  at  equispaced  frequencies  over  one  period  of  the 
desired  response 

H(k)  =  H(e>9);  9  =  2-Kk/N.  (3.39) 

If  the  number  of  samples,  N ',  is  even,  they  need  to  be  transformed  using 
the  following  relationship,  prior  to  proceeding  with  Step  2 

H{k)  =  G{k)ej2nnklN  (3.40) 

such  that 

G(0)       =  H(0) 
G(N/2)  =  0 

G(k)      =  -G(N  -  Jfe);  k  =  1, 2, . . . ,  (N/2)  -  1. 
Step  2:     Determine  the  Fourier  coefficients  of  the  desired  filter  by  finding  the 
IDFT  of  these  sample  values.    In  other  words,  determine  the  impulse 
response. 

h(n)  =  (1/N)  J]  H(k)ej2irnk/N  (3.41) 
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Step  3:    Again,  it  may  be  necessary  to  use  an  appropriate  window  function, 

h(n)  =  h(n)  •  w(n).  (3.42) 

Step  4:    The  Fourier  coefficients  thus  determined  correspond  to  the  weighted  filter 

impulse  response,  h(n),  and  the  filter  is  realized  nonrecursively  by  the 

difference  equation 

N-i 

y(n)  =  Yl  b^n  ~  fc)  (3-43) 


fc=0 

where 

bk  =  h(n). 

Step  5:    Verification  of  the  filter  design  is  accomplished  by  determining  the  inter- 
polated frequency  response,  H(e^e),  resulting  from  the  use  of  the  above 
filter  coefficients.    This  frequency  response  is  compared  to  the  original 
desired  frequency  response,  Hd(^9),  to  see  if  it  is  a  reasonable  approxi- 
mation. 
This  procedure  is  illustrated  by  the  design  of  a  bandpass  filter.  First,  an  odd 
number  of  frequency  samples  will  be  used  and  then  an  even  number. 
Example  3.1      Frequency  Sampling  for  N  Odd 

A  bandpass  filter  is  desired  with  the  ideal  frequency  response  characteristic 
illustrated  below.  The  number  of  frequency  samples  to  be  used  is  N  =  255. 
Sampling  gives 

9k  =  2rrk/N 


where  N  =  255 

(3.44) 
k  =  0,1,...,  254 


or   A0  =  2;r/iV  =  0.0246. 
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(e       ) 


1.571       2.617         R  3.666       4.712 


2a-      6 


Figure  3.13.      Desired  Bandpass  Filter  Frequency  Response 

Thus, 

(-  1.0;  k  =  64, . . . ,  106     and  it  =  140, . . . ,  101 

JJ(fc)  =  I  0.0;  ifc  =  0, . . . ,  63         and  it  =  107, . . . ,  148 
I  and  fc=  192,..., 254  . 

Taking  the  ID  FT  of  the  above  frequency  sample  values  and  truncating  h(n)  to  51 

terms  with  a  rectangular  window,  produces  the  interpolated  frequency  response  of 

Figure  3.14a,  while  Figure  3.14b  illustrates  the  response  obtained  using  a  Hamming 

window. 

Example  3.2      Frequency  Sampling  for  N  Even 

Repeat  the  previous  example,  with  N  =  256  frequency  samples. 

Sampling  gives 


where 


0k  =  2zk/X 

N  =  256 
A:  =  0,1,... 


A0  =  2»/.V  =  0.0245. 
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BANDPASS  FILTER  (N  =  255) 
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Figure  3.14a.     Unwindowed  Bandpass  Filter  Frequency  Response 
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BANDPASS  FILTER  (N  =  255 


1              ! 

i              1         / 

\! 

!             i        / 

\ 

/ 

\    ' 

/ 

• 

!       / 

i              1       / 

\ 

!                   / 

!          1 

:         / 

\ 

!                 1      / 

\ 
\ 

/ 

\ 

0.00000    0.62832      1.25664       1.88496 

FREQUENCY,  THETA 


2.51321 


3.14160 


HAMMING  WINDOW 


Figure  3.14b.     Windowed  Bandpass  Filter  Frequency  Response 
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Thus, 


f  1.0; 
=  <  0.0; 


1.0;     fc  =  64,...,106     and  fc  =  150,...,  192 

{  0.0;      k  =  0, . . . ,  63        and  k  =  107, . . . ,  149 

and  k  =  193, . . . ,  255 


But  recall  for  N  even,  the  sample  values  H(k),  must  be  transformed  to  elimi- 
nate the  unwanted  imaginary  component.  Applying  Equation  (3.40)  gives 

H(k)  =  G(k)e^k/256 


such  that, 


For  example 


G(0)  =  H(0)  =  0 

G(256/2)  =  G(128)  =  0 

G(k)  =  -G(256  -  Jfe);  k  =  1, 2, . . . ,  127. 

G(65)  =  H(65)  =  1.0  =  -G(191)  =>  G(191)  =  -1.0 


and  the  transformed  #(65)  and  if  (191)  are 

#(65)  --=  i.Oe'2^65)/256 

tf(191)  =  -1.0ej2,r(191)/256. 

Using  the  transformed  H(k)  values,  the  IDFT  is  determined  to  find  the  impulse 
response,  which  is  again  truncated  to  51  terms.  Figures  3.15a  and  3.15b  show  the 
unwindowed  and  windowed  frequency  responses,  respectively.  It  should  be  noted 
that  when  using  a  Hamming  window  for  an  even  number  of  terms  (in  this  case  52) 
the  equation  for  the  window  function  (Equation  3.19)  should  be  modified  as  follows 


w(n)  =  0.54 +  0.46  cos 
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Figure  3.15a.     Unwindowed  Bandpass  Filter  Frequency  Response 
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BANDPASS  FILTER  (N  =  256) 
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Figure  3.15b.     Windowed  Bandpass  Filter  Frequency  Response 
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G.    TRANSITION  POINTS 

As  discussed  in  the  previous  section,  frequency  sampling  involves  taking  N 
equispaced  samples  of  a  desired  filter  frequency  response,  Hd(^9)-  The  IDFT 
algorithm  is  used  to  approximate  the  unit  sample  response,  /i(n),  which,  for  the 
case  of  nonrecursive  filters,  is  equivalent  to  the  Fourier  series  coefficients.  The 
coefficients,  in  turn,  are  used  to  design  a  filter  that  approximates  the  original 
desired  continuous  frequency  response. 

A  filter  designed  in  this  manner  has  a  frequency  response  that  is  equal  to 
the  desired  frequency  response  at  the  frequencies  of  the  sampled  values,  however, 
the  response  can  deviate  significantly  from  the  desired  response  at  frequencies  in 
between  these  sample  values. 

A  method  that  can  be  used  to  smooth  the  frequency  response  involves  the 
use  of  transition  samples  between  the  passband  and  stopband  of  the  desired  fre- 
quency response  [9].  The  values  of  these  transition  samples  are  selected  based  on 
minimization  of  the  ripple  in  the  passband  and/or  stopband,  or  minimization  of 
the  maximum  sidelobe  of  the  frequency  response.  In  other  words,  a  minimization 
algorithm  is  used  to  find  the  optimum  values  for  the  transition  samples  based  on 
the  selected  minimization  criterion. 

As  before,  an  example  will  be  given  to  illustrate  this  technique,  but  first  a 
summary  of  the  steps  involved: 

Step  1:  Sample  the  desired  continuous  frequency  response  at  N  equispaced  fre- 
quencies where  N  is  the  number  of  Fourier  coefficients  that  will  be  used 
to  approximate  the  desired  filter  response.  The  spacing  between  the 
frequencies  is  A0  =  2it/N. 
Step  2:  Determine  the  impulse  response  h(n)  (Fourier  coefficients)  by  finding  the 
Inverse  Discrete  Fourier  Transform,  IDFT,  of  the  sample  values. 
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Step  3:    Using  the  coefficient  values  determined  in  the  previous  step,  find  the  con- 
tinuous frequency  response  and  compare  with  the  desired  filter  response. 
Step  4:    Use  the  minimization  algorithm  to  adjust  the  transition  coefficient  values, 
in  order  to  obtain  as  close  a  match  as  possible  between  the  desired  filter 
response  and  that  obtained  through  the  design  procedure. 
The  minimization  algorithm  as  used  in  [9],  is  based  on  minimizing  the  maxi- 
mum sidelobe  of  the  frequency  response.  Tables  were  generated  for  lowpass  filters, 
bandpass  filters  and  wide-band  differentiators  of  varying  bandwidths;  1,  2,  or  3 
transition  points;  and  odd  and  even  values  of  the  number  of  frequency  samples,  N. 
Table  3.1  is  a  reprint  of  subsets  of  these  tables  (due  to  Reference  9)  for  lowpass 
filter  design,  using  1,  2,  or  3  transition  points.  The  number  of  frequency  samples  is 
N  =  15,  the  column  labeled  minimax  refers  to  the  maximum  sidelobe,  and  Type-1 
data  means  that  the  first  frequency  sample  is  taken  at  0  =  0. 
Design  Example 

In  this  example  the  goal  is  to  design  a  lowpass  filter  with  three  frequency 
samples  in  the  passband  symmetric  about  the  origin  (BW  =  3),  and  a  total  of  15 
frequency  samples  (JV  =  15).  The  effects  of  using  0,  1,  2,  and  3  transition  points 
will  be  investigated. 

The  values  of  the  frequency  samples,  transition  point  values,  and  impulse 
responses,  h(n),  are  summarized  below  for  all  four  cases;  the  transition  point  values 
were  obtained  from  Table  3.1. 
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TABLE  3.1 

LOWPASS  FILTER  DESIGN 

(TYPE-1  DATA,  N  ODD) 

ONE  TRANSITION  COEFFICIENT 

(due  to  reference  [9]) 
BW        Minimax  7\ 


TV  =  15 

1 

-42.30932283 

0.43378296 

2 

-41.26299286 

0.41793823 

3 

-41.25333786 

0.41047363 

4 

-41.94907713 

0.40405884 

5 

-44.37124538 

0.39268189 

-56.01416588     0.35766525 


TWO  TRANSITION  COEFFICIENTS 

BW         Minimax  Ti  T2 


JV 

=  15 

1 

-70.60540585 

0.09500122 

0.58995418 

2 

-69.26168156 

0.10319824 

0.59347118 

3 

-69.91973495 

0.10083618 

0.58594327 

4 

-75.51172256 

0.08407593 

0.55715312 

5 

-103.46078300 

0.05180206 

0.49917424 

THREE  TRANSITION  COEFFICIENTS 

BW         Minimax  Tx  T2  T3 


JV  =  15 

1  -94.61166191  0.01455078 

2  -104.99813080  0.01000977 

3  -114.90719318  0.00873413 

4  -157.29257584  0.00378799 


0.18457882  0.66897613 

0.17360713  0.65951526 

0.16397310  0.64711264 

0.12393963  0.60181154 
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Case  1  -  0  Transition  Points 


H(k) 


N    =    15 


-i    _6    _5    _4    -3    _2    -i 


12       3       4       5       6       7 


Figure  3.16a.  Case  1  -  Frequency  Samples 

k     ->     0  ±1  ±2  ±3     ±4  ±5  ±6  ±7 

fl"(fe)     ->     1  1  1  0        0  0  0  0 

h{n)     -♦  0.333  0.278  0.142  0.0  -0.077  -0.067  0.0  0.058 
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Case  2  -  1  Transition  Point 
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T, =    0.410474 


12       3       4       5       6       7  k 


Figure  3.16b.      Case  2  -  Frequency  Samples 


k  -> 
H{k)  -, 
h(n)     -^ 


0  ±1         ±2         ±3  ±4  ±5  ±6         ±7 

1  1  1  0.410        0  0  0  0 
0.388     0.295     0.098     -0.044     -0.061     -0.012     0.017     0.014 
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Case  3  -  2  Transition  Points 


H(k) 


N    =    15 

T\=    0.100836 


■1    -6    -5    -4 


Figure  3.16c.      Case  3  -  Frequency  Samples 

k     —     0  ±1         ±2         ±3  ±4  ±5         ±6  ±7 

H{k)     -*     1  1  1  0.586        0.101        0  0  0 

h(n)     ->     0.425     0.300     0.066     -0.059     -0.041     0.005     0.013  0.004 
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Case  4  -  3  Transition  Points 


H(k) 


N    = 

15 

T,  = 

0.008734 

T2  = 

0.163973 

.T3 

T3  = 

0.647113 

3       4       5       6       7  k 


k     ->     0  ±1         ±2         ±3  ±4  ±5         ±6  ±7 

H(k)     ->     1  1  1  0.647        0.164        0.009     0  0 

h{n)     -*     0.443     0.301     0.050     -0.062     -0.032     0.008     0.010  0.002 

Figure  3.16d.      Case  4  -  Frequency  Samples 
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Figure  3.17  depicts  the  results  of  frequency  sampling  design  using  transi- 
tion points.  As  expected,  the  amplitude  of  the  passband  and  stopband  ripple  is 
reduced  with  an  increasing  number  of  transition  points;  however,  the  tradeoff  is 
that  the  sharpness  of  the  cutoff  is  reduced.  This  example  illustrates  the  usefulness 
of  minimization  algorithms  in  smoothing  interpolated  filter  frequency  responses 
obtained  with  the  frequency  sampling  technique.  As  stated  earlier,  Reference  9 
contains  tables  for  lowpass  and  bandpass  filter  design,  as  well  as  wideband  differ- 
entiators. While  quite  extensive,  a  designer  may  want  parameters  that  have  not 
been  tabulated.  In  this  case,  approximate  values  of  the  transition  coefficients  can 
be  derived  from  the  tabulated  values  given,  using  linear  interpolation. 

H.    DESIGN  OF  AN  INTEGRATOR 

In  Section  A  analytical  methods  were  used  to  design  a  differentiator.  The  use  of 
these  methods  to  design  an  integrator,  however,  does  not  yield  an  easily  obtainable 
solution.  As  will  be  shown  shortly,  frequency  sampling  solves  this  dilemma  by 
providing  a  straightforward  design  technique  that  produces  extremely  good  results. 

Bandpass  Integrator  Design  Example 

A  bandpass  integrator  with  the  following  frequency  response  characteristic 
is  desired.  Both  of  the  aforementioned  design  methods  will  be  applied  to  the  design 
problem  to  illustrate  the  advantage  of  the  frequency  sampling  technique. 

a.    Method  1  -  Analytical 

In  the  5-domain  an  integrator  is  described  by  the  transfer  function 
H(s)  =  1/s.  The  frequency  response  is  thus:  H(jw)  =  1/jw  =  {jw)~l .  Converting 
the  frequency  response  to  digital  frequencies 

H(jw)  =  1/jw  =>  H(ej0)  =  -jT/0-  T  =   period.  (3.46) 


124 


-■ 

1 

! 

o 

/?-\ 

LEGEND 

0  TRANS  PT 

1  TRANS  PT 

2  TRANS  PT 

3  TRANS  PT 

— 

05 

^ — "                 V\N           i 

ep 

!      \\   i 

d~ 

!       \  V  \ ! 

Q  d~ 

!        \  ■  \( 

t*  "R 

tU  d~ 
Z  « 

\\f\           i             1 

I         \    i  \\         I             1 

"*5    O 

\  i-  \\        ! 

©  ~ 

\  !  ■   V\       ! 

d  " 

© 

d 

1              \ 

/       \/^^^\ 

/\/ 

0..00000        0.62832  1.25664  1.88-196 

FREQUENCY,  THETA 


:.5132< 


3.14160 


Figure  3.17.     Frequency  Responses  for  Varying  Numbers 
of  Transition  Points 
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K  -8, 


Hn(e       ) 


Figure  3.18.     Ideal  Integrator  Frequency  Response  Characteristic 

To  determine  the  impulse  response,  h(n), 

fc(n)  =  (1/2tt)  /        tfD  (e^)  e^-c/0  4-  1/2t  /     tfD  (eje)  e>9"dO 
J-97  Jex 

"  dd 


& 


h(n)  =  (l/2rr)  f       (-JT/8)  eJe"d9  +  1/2*  f     l/6ej 
Since  :   f{l/x)eax 
h(n)  =  -;T/2;r 


Jx  =  lnX  +  l!+2^  +  ^l!+- 


,  ,  jnfl  ,   (jn)W   ,   (jn)303    , 
ln  9  +  —  +  T^T  +  T"3T  + 


A  considerable  number  of  steps  later  yields 


jn0   ,    On)2^2    |    (jn)3^    | 


t      +      9.0) 


3-3! 


h(n)  =  T/x 


n(92-el)  +  —  (el-el)+...  +  —  (8?-e?) 


(3.47) 


where  n  is  the  number  of  coefficients  desired  for  the  implementation  of  the  bandpass 
integrator. 

The  analytical  expression  for  the  impulse  response  of  the  bandpass 
integrator,  Equation  (3.40),  is  cumbersome,  especially  when  the  integrator  design 
is  of  a  high  order.  Therefore,  an  alternative  method  is  desirable. 


126 


b.    Method  2  -  Frequency  Sampling 

Recall  that  the  expression  for  the  desired  frequency  response  is 

HD  (e*)  =  -jT/0.  (3.48) 

Letting  T  =  1,  one  hundred  and  one  frequency  samples  will  be  taken  of  the  de- 
sired bandpass  integrator  with  a  bandwidth  of  OAtt  to  1.67T.  In  other  words,  101 
frequency  samples,  H(k),  of  the  magnitude  and  phase  characteristic  (Figure  3.19) 
will  be  taken  using  the  following  increment  of  the  digital  frequency  0 

A0  =  2wk/N  =  0.062A;  for  TV  =  101 
H{k)  =  H  (ej°)  for  k  =  -50, . . . ,  +50. 

0=O.O62Jfc 

The  101-coefficient  impulse  response  is  obtained  by  determining  the 
IDFT  of  these  frequency  sample  values,  and  a  21-coemcient  rectangular  window  is 
applied  yielding  the  following  bandpass  integrator  coefficients. 


R 

h{n) 

R 

h(n) 

0 

0.0974 

11 

0.4813 

1 

0.0793 

12 

0.2391 

2 

0.0946 

13 

0.2248 

3 

0.0439 

14 

0.0818 

4 

0.0362 

15 

0.0562 

5 

-0.0562 

16 

-0.0362 

6 

-0.0818 

17 

-0.0439 

7 

-0.2248 

18 

-0.0946 

8 

-0.2391 

19 

-0.0793 

9 

-0.4813 

20 

-0.0974 

10 

0.000 

Finally,  the  interpolated  frequency  response  is  obtained.  As  can  be 
seen  in  Figure  3.20,  the  results  are  good.  (Solid  dots  indicate  the  ideal  frequency 
response.)  Also  included  in  Figure  3.21  are  the  results  obtained  when  a  Hamming 
window  is  used,  while  Figures  3.22  and  3.23  are  the  results  of  using  41  coefficients. 
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Figure  3.19.      Ideal  Bandpass  Integrator  Frequency  Response  Characteristic 
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This  example  further  demonstrates  the  versatility  and  usefulness  of 
the  frequency  sampling  technique.  It  allows  for  the  design  of  filters  that  may  be 
unattainable  using  traditional  analytic  methods,  and  is  thus  an  indispensable  tool. 
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Figure  3.20.     Unwindowed  21-Coefficient  Bandpass  Integrators 
(Solid  dots  indicate  ideal  frequency  response) 
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Figure  3.21.     Windowed  21-Coefficient  Bandpass  Integrator 
(Solid  dots  indicate  ideal  frequency  response) 
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Figure  3.22.     Unwindowed  41-Coefficient  Bandpass  Integrator 
(Solid  dots  indicate  ideal  frequency  response) 
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Figure  3.23.     Windowed  41-Coefficient  Bandpass  Integrator 
(Solid  dots  indicate  ideal  frequency  response) 
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IV.    COMPUTER-AIDED  DESIGN 

The  basis  of  all  computer-aided  design  (CAD)  techniques  is  optimization.  This 
is  to  say,  a  desired  filter  frequency  response  is  approximated  by  a  particular  filter 
whose  coefficients  are  to  be  determined.  The  accuracy  of  the  approximation  is 
evaluated  according  to  some  criterion,  usually  an  error  function,  that  indicates 
how  large  a  disparity  exists  between  the  desired  filter  frequency  response  and  the 
approximated  filter  frequency  response.  Variable  parameters  of  the  approximating 
function  are  then  "adjusted"  to  optimize  the  filter  design  in  terms  of  this  criterion. 

A.    REMEZ  EXCHANGE  ALGORITHM 

An  extensively  used  computer-aided  design  technique  for  linear  phase  FIR 
filters  is  the  Remez  exchange  algorithm.  The  mathematical  basis  for  this  algorithm 
is  the  weighted  Chebyshev  approximation. 

A  summary  of  the  approximating  and  error  functions  for  this  algorithm  follows. 
It  has  been  shown  in  [1]  and  [13],  that  the  frequency  responses  for  the  four  cases 
of  linear  phase  filters  -  i.e.,  even  or  odd  symmetry  with  an  even  or  odd  number  of 
terms  -  can  be  written  in  the  form: 

H  (e?9)  =  t-W-WjWWH  (e>'«)  (4.1) 

where  H  (ejtf)  is  a  real-valued  function  used  to  approximate  the  desired  niter's 
magnitude  specifications,  and  the  remaining  terms  approximate  the  desired  phase. 
Table  4.1  gives  the  values  L,  and  the  form  of  H  (e7'*)  for  all  four  cases. 


134 


TABLE  4.1 
FREQUENCY  RESPONSES  FOR  LINEAR  PHASE  FILTERS 

(due  to  reference  13) 

L    H  (eg) 
Case  1  -  N  odd 

(N-l)/2 

Symmetrical  0  ^      a(n)  cos(n#) 

n=0 

Impulse  Response 


Case  2  -  JV  even 

iV/2 

Symmetrical  0      £  6(rc)  cos  [9(n  -  1/2)] 


n=l 


Impulse  Response 


Case  3  -  JV  odd 

(N-l)/2 

Anti-symmetric  1  ^      c(n)sin(n#) 

n=l 

Impulse  Response 


Case  4  -  iV  even 

AT/2 

Anti-symmetric         1       £]  cf(n)  sin[#(n  —  1/2)] 

n=l 

Impulse  Response 

Using  trigonometric  identities,  the  expressions  for  H  (e^e)  in  Table  4.1  can  be 
rewritten  in  the  following  form: 

H  (e")  =  Q  (e?9)  P  (e*)  (4.2) 

where  Q  (e7*)  is  a  fixed  function  of  frequency,  0,  and  P  (e-7^)  consists  of  a  sum  of 
weighted  cosine  terms  (see  Table  4.2). 


135 


TABLE  4.2 
FREQUENCY  RESPONSES  FOR  LINEAR  PHASE  FILTERS 

(due  to  reference  13) 

0(e")  P(e>«) 


(N-l)/2_ 
Case  1  1  ^      <Kn)  cos(n#) 

n=0 
(N/2)-l„ 
Case  2     cos(0/2)  X)      6(n)cos(n0) 

n=0 
(N-l)/2 

Case  3       sin(#)  ^      c(n)cos(n#) 

n=0 

(N/2)-l    _ 

Case  4     sin(0/2)  £      d(n)cos(n#) 

n=0 


where 


a{n>-\a(n)  =  2h(^-n: 


for  n  =  1,2,.       is=1 


2 


r  6(1)  =  6(0)  +  1/26(1) 
b(n)=  I  6(n)  =  1/2  [&(n  -  1)  +  6(n)l     for  n  =  2,3,... ,  f  -  1 

l6(iV/2)  =  l/26(f -1) 


c(l)  =  c(0)-l/2c(2) 

c(n)  =  1/2  [c(n  -  1)  -  c(n  +  1)]     for  n  =  2,3, . . . ,  ^f1- 
c(n)=<j  c(ivfi_l)  =  1/2g(A^i_2) 

c(^i)  =  l/2c(^i-l) 


d(l)  =  J(0)  -  1/2J(1) 
^  <f(n)  =  1/2  [J(n  -  1)  -  J(n)]      for  n  =  2,3, . . . ,  f  -  1 
<f(f)  =  l/2J(f-l) 
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Since  Q  (e-*6)  is  a  function  of  frequency  only,  it  can  be  seen  that  the  approx- 
imating function,  H  (eJ°),  can  only  be  optimized  in  terms  of  P  (e-7*);  that  is,  the 
dependency  of  H  (e-7^)  on  the  filter  coefficients  is  contained  in  P  (e-7  ).  Thus,  the 
coefficients  of  P  (e^9)  can  be  varied  to  achieve  an  optimum  filter  design,  and,  the 
true  approximating  function  can  be  generalized  for  all  cases  as, 


P{e?9)  =  ^a(n)cos(n0) 


(4.3) 


n=0 


where  L  =  (N  —  l)/2  or  (N/2)  —  1  depending  on  which  case  is  being  considered, 
and  the  a(n)  are  the  weights  from  which  the  filter  coefficients  can  be  determined. 
As  stated  in  the  introduction  to  this  chapter,  a  measure  of  how  well  the  de- 
signed filter  frequency  response  approximates  the  desired  filter  frequency  response 
is  required.  The  weighted  Chebyshev  approximation  uses  an  error  function  defined 
as  follows: 

E($)  =  W{8)  ]HD  (e?e)  -  H  (e>*)]  (4.4) 

where 

Hd  (eJ  )  =  the  desired  frequency  response 

H  (eJ  )  =  the  designed  frequency  response 

W(0)  =  weighting  factor 

E{9)  =  error 
In  order  to  see  the  relationship  between  the  filter  coefficients  and  the  error, 
Equation  (4.4)  can  be  rewritten  in  terms  of  the  functions  Q  (eJ*)  and  P  (ejd). 

E(9)  =  W(9)  [HD  {e?9)  -  Q  (e")  P  (e*)] 
HD(e^) 


W{9)Q  (e") 


(e") 


-p(*n 


Letting 


W{9)  =  W(9)Q  (e^)  and  HD  (e>8)  =  HD  (e*#)  /Q  (<J 


(4.5) 


(4.6) 
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yields 

E(9)  =  W(9)  [HD  (e*)  -  P  (e?6)]  (4.7) 

Thus,  the  Chebyshev  approximation  that  is  performed  using  the  Remez  exchange 
algorithm  can  be  stated  as  follows: 

Find  the  set  of  filter  coefficients  (determined  from  the  values  of  a(n))  that 
minimizes  the  maximum  absolute  value  of  the  error,  E(9),  over  the  frequency  range 
of  interest. 

n 

e  \E(6)\ 


(4.8) 


-E optimum  —  min 

(coefficients) 

At  this  point,  a  discussion  of  the  weighting  function,  W{9),  is  in  order.  The 
purpose  of  this  weighting  function  is  to  ensure  a  small  tolerance  for  error  in  critical 
frequency  ranges. 

If  W{9)  is  large,  this  means  a  large  deviation  between  the  desired  frequency 
response,  Ho  (e^*))  and  the  designed  frequency  response,  H  (e-7*),  cannot  be  toler- 
ated. Looking  at  Equation  (4.4)  we  see  that  if  W{9)  is  large,  the  difference  between 
the  desired  and  designed  frequency  responses,  \Hd  (e-7*)  —  H  (eJ(9)  has  to  be  small 
to  keep  the  weighted  error  small. 

Conversely,  if  W(9)  is  small,  the  difference  Ed  (^e)  —  H  (e-7*)  can  be  larger 
and  still  meet  the  error  criterion.  Small  values  for  W{9)  would  be  used  in  frequency 
bands  where  close  approximation  to  the  desired  frequency  response  is  not  critical. 

A  copy  of  a  program  employing  the  Remez  exchange  algorithm  [16]  lias  been 
included,  and  its  use  is  best  illustrated  with  an  example.  Before  proceeding  with 
the  example,  the  salient  features  of  the  FIR  Linear  Phase  Filter  Design  program 
can  be  summarized  as  follows. 
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1.    Main  Program  Functions 

Manage  the  Input  - 

NFILT  =  filter  length  in  terms  of  samples,  3  <  NFILT  <  NFMAX 
NFMAX  =  128,  but  can  made  greater  or  smaller  by  the  user 
JTYPE  =  filter  type 

1  =  Multi-Passband/Stopband 

2  =  Differentiator 

3  =  Hilbert  Transformer 

EDGE  =  number  of  frequency  bands,  specified  by  0iower  and  0upper 

(  maximum  allowed  is  10  ) 
FX  =  desired  frequency  response  magnitude  in  each  band,  \Hd  {eJ&)  \ 
WTX  =  positive  weight  function,  W(9),  in  each  band 
LGRID  =  Grid  Density  (the  grid  density,  default  value  is  16) 

Set-up  appropriate  approximation  problem,  based  on  the  desired  frequency 
response,  Hd  (ejd)»  and  the  weights,  W{9). 

Yield  Output  - 

o  The  coefficients  of  the  best  impulse  response  obtained  from  the  best  cosine 
approximation. 

rmax 


\E{9)\ 


o  The  optimal  error  I  min 

max 

o  The  extremal  frequencies  where  E{9)  =  9    \E(9)\ 
2.    Most  Important  Subroutines 

EFF  -  defines  the  desired  frequency  response,  Hd  (ejd) 
WATE  -  defines  the  weight  function,  W{9) 
REMEZ  -  calculates  the  best  Chebyshev  approximation  of  the  desired  frequency 
response 


139 


3.  Important  Note 

It  should  be  noted  that  the  program  output  is  the  impulse  response,  h(n). 
The  user  must  supply  his/her  own  program  to  calculate  and  plot  the  filter  frequency 
response,  which  is  necessary  to  determine  the  suitability  of  the  filter  design. 

4.  Input  Format 

NFILT      JTYPE      #  Bands      LGRID 
EDGE  (Band  Edges) 
FX  (Magnitude  in  Bands) 
WTX  (Weights  in  Bands) 

5.  Bandpass  Filter  Design  Example 

We  wish  to  design  a  bandpass  filter  that  meets  the  following  specifications: 
Passband  with  a  gain  of  one  for  the  frequency  interval  10  kHz  to  15  kHz, 
with  a  system  sampling  frequency  of  50  kHz. 

•  The  digital  frequency  band  is: 

2tt(104) 

2t(1.5)  (104) 
'•=        5(10<)        =0-5" 

•  Normalizing  so  that  0  =  tt  corresponds  to  a  normalized  frequency  of  0.5. 

0.47r  _   6j_ 
w     ~  0.5 

o.6tt     eu 


.  0e  =  0.2 

7T  0.5 


0.5  ^  =  °-3 
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•  Designate  Band  Edges,  Magnitudes,  Weights:  (the  weights  were  arbitrar- 
ily made  equal  for  this  example) 

Band  1 

Edges  -  0.0  to  0.19 

Magnitude  -  0.0 

Weight  -  10.0 

Band  2 

Edges  -  0.2  to  0.3 
Magnitude  -  1.0 
Weight  -  10.0 

Band  3 

Edges  -  0.31  to  0.5 
Magnitude  -  0.0 
Weight  -  10.0 

•  Select  Filter  Length  (  number  of  coefficients  ) 

iV  =  21 
6.    Sample  Input  File 
00021     00001     00003     0 
0.0        0.19      0.2       0.3  Band  Edges 

0.31      0.5 

0.0        1.0        0.0  Magnitudes 

10.0      10.0      10.0  Weights 
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7.    Sample  Output  File 

The  output  file  below  was  generated  using  the  above  input  file.  Njte 
the  output  consists  of  the  filter  impulse  response,  the  input  values  used  and  the 
deviation  between  the  desired  and  approximated  frequency  responses. 


pj'ITf    T[*T,rTL^c'    E ESPOUSE     (ftt| 
LINEAP' "HAEE    OIOTm\L  '  ^TT.T^0'  6*STGN 
°EME7    ErCFVIG17    ftfiO^ITHI 

P1NDnASS    ETLTEP 

"TITE75    LENGTH    =       71 

*****      T1"»fJT  f^'?      orqr>QMq^      ***** 
r!(      1)       =      "^.11^^0311X^  +  00      =      q   (      9  1  ) 

Tf  f  2)  =  -0.  1037?3^  1"-0?  =  R(  to) 
w(  3)  =  0.110971011? +00  =  '-j(  io> 
H(     4)      =     -0.23107301^-0  3     =     it  (      1fl) 

h(   s>    =  -o. 1^^307^+00   =   a  (    17^ 

»(  M  =  -0.  2'?3">17U*E-03  =  »f  16' 
H       7)     =       0.  18^3  18f^a^  +  oo    =    jj  )     1S« 

Hf    «j  =    -0.  1?8  1S3U^^-0  3  =  '■»(  14' 

h(    o\  -    -0. 2207nq21r+OQ  =  n(  n 

H    10)  =    -0.  0^6  1FO03F-0U  =  U(  12' 
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8.    Results 

Figures  4.1  to  4.3,  illustrate  the  frequency  responses  for  filter  lengths  N  = 
21,41  and  61,  respectively. 

As  expected,  the  frequency  response  improves  with  an  increasing  number 
of  coefficients.  The  cutoff  is  sharper  and  the  magnitude  more  closely  approximates 
0  dB  in  the  passband.   As  with  other  design  methods,  windows  may  also  be  used 
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to  improve  the  frequency  response  characteristic.  Figures  4.4  to  4.6,  show  the  effect 
of  changing  the  band  edges  as  follows,  for  N  =  21  coefficients. 
Figure  4.4  -  Band  1         0.0  to  0.19 

Band  2         0.2  to  0.3 

Band  3         0.31  to  0.5 

Figure  4.5  -  Band  1         0.0  to  0.15 
Band  2         0.2  to  0.3 
Band  3         0.35  to  0.5 

Figure  4.6  -  Band  1  0.0  to  0.1 
Band  2  0.2  to  0.3 
Band  3         0.4  to  0.5 

Figure  4.4  represents  the  original  cutoff  frequencies,  while  Figures  4.5  and 
4.6  show  the  results  of  not  selecting  the  stopband  frequency  cutoff  close  enough  to 
the  passband  frequency  starting  point,  i.e.,  the  passband  is  broadened  beyond  the 
desired  0.2  to  0.3  specifications. 

Thus,  when  selecting  the  stopband  cutoff  frequencies  one  should  choose 

values  as  close  to  the  passband  frequencies  as  possible.  It  is  of  importance  to  note 

that  the  program  does  not  allow  for  the  upper  stopband  frequency  equaling  the 

lower  passband  frequency.   In  other  words  the  following  parameters  would  not  be 

acceptable: 

Band  1     0.0  to  0.2 

Band  2     0.2  to  0.3 

Band  3     0.3  to  0.5 
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Figure  4.1.     21-Coefficient  Remez  Bandpass  Filter  Design 
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Figure  4.2.     41-Coefficient  Remez  Bandpass  Filter  Design 
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Figure  4.3.     61-Coefficient  Remez  Bandpass  Filter  Design 
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Figure  4.4.     21-Coefficient  Remez  Bandpass  Filter  Design 
(Band  edge  spacing  0.01) 
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(Band  edge  spacing  0.05) 
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Figure  4.6.     21-Coefficient  Remez  Bandpass  Filter  Design 
(Band  edge  spacing  0.1) 
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Figure  4.7.     21-Coefficient  Remez  Bandpass  Filter  Design 
(Band  weighting  1:1:1) 
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Figure  4.8.      21-Coefficient  Remez  Bandpass  Filter  Design 
(Band  weighting  10  :  1  :  10) 
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Figure  4.9.     21-Coefficient  Remez  Bandpass  Filter  Design 
(Band  weighting  10  :  10  :  10) 
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Finally,  Figures  4.7  to  4.9  are  the  result  of  varying  the  weights  of  the 
original  filter  specifications  as  follows  for  N  =  21. 

Figure  4.7  -  Band  1  weight  =  1.0 
Band  2  weight  =  1.0 
Band  3     weight  =  1.0 

Figure  4.8  -  Band  1     weight  =  10.0 
Band  2     weight  =  1.0 
Band  3     weight  =  10.0 

Figure  4.9  -  Band  1  weight  =  10.0 
Band  2  weight  =  10.0 
Band  3     weight  =  10.0 

Again,  as  expected,  the  frequency  bands  with  the  heavier  weight  assigned 
more  closely  approximated  the  desired  frequency  response,  Figure  4.8.  Here  it 
should  be  noted  that  the  weights  are  taken  into  account  relative  to  one  another. 
Figures  4.7  and  4.9  are  identical  because  in  both  cases  the  relationships  between 
the  weights  is  1:1:1. 

In  this  example,  we  have  seen  that  the  Remez  exchange  algorithm  provides 
an  effective  way  to  design  linear  phase  FIR  filters.  However,  a  primary  drawback 
to  this  procedure  is  that  the  CPU  time  required  grows  quite  rapidly  (on  the  order 
of  N)  with  the  filter  order  N. 

For  example,  fifteen  CPU  minutes  are  required  for  the  design  of  a  lowpass 
filter  with  N  =  512  using  a  CDC  6500  [15].  This  is  the  time  required  for  one 
iteration  of  the  procedure.   Typically,  more  than  one  iteration  is  required  to  find 
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technique  to  reduce  the  CPU  time  required  is  desirable.  Reference  15  presents  a 
method  whereby  the  pA-operties  of  a  high-order  filter  can  be  extrapolated  from  a 
lower-order  filter. 

B.    METHOD  FOR  THE  DESIGN  OF  HIGH-ORDER  LINEAR  PHASE 
FIR  FILTERS  BASED  ON  A  LOW-ORDER  PROTOTYPE 

A  detailed  explanation  of  this  method  will  not  be  presented,  however,  a  sum- 
mary of  the  general  idea  behind  this  technique  follows.  The  interested  reader  is 
directed  to  Reference  15  for  details  of  its  application.  It  should  be  noted  that  the 
term  "high-order"  refers  to  filters  with  orders  approaching  N  =  2048. 

The  underlying  basis  for  the  high-order  filter  design  technique  is  the  observa- 
tion that  in  high-order,  multi-passband/stopband  filters,  extremal  frequencies  (i.e., 
ripple  frequencies)  in  the  broad  part  of  these  bands  are  more  evenly  distributed 
than  in  the  transition  regions. 
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Figure  4.10.      Desired  High-Order  Filter 

Figure  4.10  represents  the  desired  frequency  response  of  a  high-order  multi- 
passband/stopband  filter;  and  it  is  noted  that  in  regions  A,  C,  and  E  there  is  a 
uniform  distribution  of  extremal  frequencies,  while  in  regions  B  and  D  the  distri- 
bution is  non-uniform. 
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The  design  of  this  filter  using  a  lower  order  prototype  involves  the  following 


steps: 


To  obtain  the  low-order  prototype,  the  uniform  extremal  frequency  regions,  A, 
C  and  E  are  "cut-out".  The  remaining  non-uniform  regions  B  and  D  are  then 
merged  to  form  the  lower  order  filter  whose  frequency  response  Hp  (ej8  )  is 
shown  in  Figure  4.11. 


Figure  4.11.      Low-Order  Filter  Prototype 

•  The  Remez  exchange  algorithm  is  then  used  to  obtain  the  filter  order,  N,  and 
the  coefficients  required  that  will  meet  the  specifications  of  H p  ( e^e  J .  This 
is  done  by  varying  the  order  and  the  weights,  until  the  order  and  weights  that 
yield  the  desired  passband  ripple  and  stopband  attenuation  are  found. 

•  Once  the  low-order  prototype  filter  is  obtained,  its  extremal  frequencies  are 
plugged  back  into  the  appropriate  regions  of  the  original  desired  high-order 
filter,  Figure  4.12. 

•  The  extremal  frequencies  are  used  as  initial  values,  and  a  final  run  of  the 
Remez  Exchange  Algorithm  is  then  made  to  obtain  the  filter  order  and  impulse 
response  of  the  high-order  filter  design.  Again,  the  filter  order  and  weights  are 
varied  until  an  optimum  design  is  obtained. 

In  conclusion  it  should  be  mentioned  that  in  Reference  21a  program  is  outlined 
that  automatically  designs  high-order  FIR  filters  using  this  procedure.  Approxima- 
tions of  filters  with  orders  up  to  N  =  4096  have  been  achieved  using  this  program. 
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Figure  4.12.     Low-order  filter  prototype  being  placed 

back  into  the  appropriate  regions  of  the 

desired  high-order  filter 

In  conclusion  it  should  be  mentioned  that  in  Reference  21a  program  is  outlined 
that  automatically  designs  high-order  FIR  filters  using  this  procedure.  Approxima- 
tions of  filters  with  orders  up  to  N  =  4096  have  been  achieved  using  this  program. 

This  section  has  illustrated  the  usefulness  and  relative  ease  of  application  of 
Computer- Aided-Design  (CAD)  techniques,  such  as  the  Remez  exchange  algorithm. 

Since  this  algorithm  is  only  applicable  to  linear  phase  FIR  filters,  it  is  worth- 
while at  this  point  to  diverge  to  a  discussion  of  a  popular  algorithm  for  the  design 
of  recursive  IIR  filters,  the  Minimum  p-Error  Design  Method. 
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C.    THE  MINIMUM  p-ERROR  DESIGN  METHOD 

Before  proceeding  with  a  discussion  of  this  technique,  it  should  be  mentioned 
that,  in  the  interest  of  brevity,  the  pertinent  equations  and  relationships  used 
are  presented  without  elaboration.  This  is  because  the  mathmatical  basis  for  the 
minimum  p-error  criterion  is  quite  extensive  and  beyond  the  scope  of  this  thesis. 
The  interested  reader  may  find  details  of  its  derivation  in  Reference  19. 

The  minimum  p-error  design  technique  consists  of  a  generalization  of  the  min- 
imum mean-squared  error  design  method,  wherein  the  optimum  filter  coefficients 
are  determined  by  minimizing  the  following  error  function  [18]. 

**M)  =  E  w (**)  & <A'  °k">  -  H^k)fp  (4.9) 

fc=i 

where 

w  (Ok)  =   weighting  factor 
H  (A,  Ok)  =   approximating  function  (actual  response) 
Hd  (Ok)  =   desired  frequency  response 

Ok  =  digital  frequencies  over  the  range  of  interest,  k 
A  =   vector  containing  the  k  independent  parameters 
(i.e.,  the  filter  coefficients) 

Note:     If  the  value  for  p  in  Equation  (4.9)  is  unity,  this  relationship  describes  the 
mean-square  error. 

The  filter  approximation  problem  can  thus  be  stated  as  follows:  Given  an 
amount  of  error  that  can  be  tolerated,  -&2P ,  find  the  set  of  k  parameters  (filter 
coefficients)  such  that  the  error  between  the  approximate  frequency  response  and 
the  desired  frequency  response  is  within  the  stated  tolerance. 
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Again,  looking  at  Equation  (4.9)  it  can  be  seen  that  for  large  values  of  p, 
the  approximate  frequency  response,  H(A,9k),  has  to  be  very  close  to  the  desired 
frequency  response,  Hoi^k),  to  remain  within  the  pre-selected  error  tolerance  value 
of  Eiv  •  Furthermore,  a  sufficiently  large  value  of  p  results  in  an  optimal  solution 
that  is  very  close  to  the  optimal  Chebyshev  (or  minimax)  solution  [19].  The  method 
of  solution  of  this  approximation  problem,  as  will  be  shown,  is  dependent  on  the 
error  tolerance  value  selected,  E2P,  the  form  of  the  transfer  function,  H(z),  and 
the  parameter  vector,  A.  First  a  discussion  of  the  transfer  function  is  presented 
because  its  form  has  a  direct  impact  on  several  important  features  in  digital  filter 
design  [18]. 

As  is  well  known,  three  possible  forms  of  the  filter  transfer  function  are  direct, 
cascade  and  parallel  realizations.  With  respect  to  errors  directly  related  to  the 
physical  structure  of  the  filter,  i.e.,  quantization  effects  and  finite  coefficient  size, 
the  use  of  the  direct  form  of  the  transfer  function  in  filter  realization  is  undesirable. 
This  leaves  the  parallel  and  cascade  forms,  with  the  cascade  form  being  selected 
because  the  zeros  of  the  transfer  function  remain  unaltered  in  the  process  of  cas- 
cading, resulting  in  well  defined  stopbands.  Conversion  of  a  filter  transfer  function 
to  the  parallel  form,  on  the  other  hand,  involves  partial  fraction  expansion,  which 
results  in  the  zeros  being  less  well  defined.  Theoretically  this  should  not  be  the 
case,  but,  when  the  parallel  form  is  implemented  digitally,  the  zeros  are  slightly 
different  than  those  of  the  original  direct  form  transfer  function  due  to  the  finite 
precision  of  the  computer. 

Using  the  cascade  form  of  the  filter  transfer  function  has  other  advantages 
including: 

•  stability  tests  of  the  filter  can  be  accomplished  without  lengthy  calculations, 
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•  the  frequency  response  is  of  a  simple  functional  form  that  readily  lends  itself 
to  insights  as  to  how  the  poles  and  zeros  (hence,  the  filter  coefficients)  impact 
the  error  function. 

Thus,  the  first  step  in  the  use  of  this  design  procedure  involves  decomposing 
the  proposed  approximating  transfer  function,  H(z),  into  cascaded  second-order 
sections. 

»=i  t=i 

where 

k0  =    a  positive  normalizing  constant,  and 

N  =   filter  order  /2. 
The  parameter  vector,  A,  therefore,  consists  of  the  multipliers  used  in  the 
cascade  structure. 

A  =  [ai1,a12,&ii,&i2,...,aif-,a2i,&i,-,&2i,-..,&o]  (4-H) 

The  following  expressions  for  the  frequency  response  and  group  delay  of  the 
filter  incorporate  the  cascade  structure  and  parameter  vector. 

•  Frequency  Response 

ntx    m_  MJ°\-  h    TT(1  +  a2,)cos^  +  qlt+i(l-Q2,)sin^ 
«{±e)-H(e>)-koil  (4.12) 


Group  Delay 


N 


T(A,8)  =  -±H(e>°)  =  Y, 


2cos0  +  &i,-  +  J2sin0 
(1  +  b2i)cos0  +  bu  +  i(l  -  b2l)  sin 


/  2cos^  +  ai, +j2sin^  N 

\(1  -1-  a2i) cos 9  +  an  +  j(l  —  a2,)sin^ 


(4.13) 


Equations  (4.12)  and  (4.13)  describe  the  frequency  response  and  group  delay 
of  the  approximating  function,  H(A,0k).    To  determine  if  the  coefficients  of  the 
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A  vector  yield  an  optimum  filter  design  based  on  the  desired  frequency  response 
Hoi^k)  and  the  error  that  will  be  tolerated  E2P(A)  involves  finding  the  partial 
derivatives  of  the  frequency  response  and  group  delay  portions  of  the  error  function, 
in  terms  of  the  filter  coefficients,  and  minimizing  them. 

In  other  words,  the  filter  coefficients  comprising  the  A  vector  are  to  be  selected 
so  as  to  minimize  the  following  partial  derivatives: 

^%^  =  t  P-.  («»)  P-  (« (A>  *•)  -  H°  (^))2'-1       (4-14) 

dau  £rj  oau 

*^*2  =  £>,(«,)  ^  (r  (A.  «*)"    HD(Sk)f>-'  (4.15) 

and  similarly  for  dE2Pa(A)/dbn  and  dE2PT(A)/dbn,  etc. 

Upon  examining  Equations  (4.12)  and  (4.13),  it  is  that  apparent  these  par- 
tial derivatives  are  not  easily  attainable  because  complex  expressions  are  involved. 
To  remedy  this  the  polar  coordinates  of  the  poles  and  zeros  are  used  for  param- 
eters, rather  than  the  original  rectangular  form.  Thus,  the  parameter  vector  and 
expressions  for  the  amplitude  and  group  delay  are  modified  as  follows, 

•  A  Parameter 

A  =  [r0i ,  0Oi ,  rpllBpU . . . ,  rtf,-,  0oi,  rpi,  9pi, . . . ,  k0]  (4.16) 

•  Frequency  Response 

°  l=\  {1  -  2rpi  cos(0  -  0~)  +  rpi2  }1/2  {l  -  2rpicos(9  +  ~$~)  +  rpi,  }1/2 

(4.17) 
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•  Group  Delay 


r(A,0)  =  £ 


l-rpicos(0-0pt) 


1  -rp,  cos  (0  +  Opt 


_  1  -  2rpi  cos  (6  -9pi)  +  rpi2       1  -  2rp,  cos  (9  +  9pt )  +  r; 


1  —  r0i  cos  (0  — 


1  -roicos(0  +  9oi 


1  -  2rot  cos  (9  -  0OI )  +  rot2        1  -  2r01  cos  (0  +  9oi )  +  r0,2 
resulting  in  the  following  partial  derivatives  of  the  error  function: 
k 


dE2pr(A] 


dr 


=  Y^P^riOk)  jr-  (r(A,Ok)  -  HD{9k)f 


droi  ^-^>  dTc 

and  for  dE2pa{A)/d9oi  and  dE2pT(A)/d90i,  etc.  where 

roi  -  cos  (9  -  90i) 


da 
drnt 


da 


+ 


1  -  2roi  cos  (9  -  9oi)  +  roj2 
roi  -  cos(9  +  9oi) 


l-2roicos(9  +  9oi 
rois'm(9  -  9oi) 


l-2roicos{9-9oi)  +  rot, 
roism(9  +  9oi) 


l-2rotcos(9  +  9ol)  +  rot2 
dr      _    (l  +  rpt-2)cos(fl-flpi)-2rpt 
&X        (l-2rpicos(0-0pt)  +  rp,2)2 


dr 
d9ot 


(l+rp,2)cos(fl  +  flpt)-2rp, 
(l-2rpicos{9  +  9pi)  +  rpi2)2 
,i(l-rpi2)sm(0-epi) 


(4.18) 


(4.19) 


(4.20) 


(l-2rp,cos(0-0pi)  +  rptO 

rpi  (l  -rp,2)sin(fl  +  flp,) 
(l-2rpicos(9  +  9pt)  +  rpt2)2' 

The  partial  derivatives  da/drpi,da/d9pi,dr/dr0i,  and  dr/d90i,  are  the  same  as 
the  above  but  with  changed  signs. 

Thus,  the  approximation  problem  amounts  to  the  minimization  of  nonlinear 
functions  (Equations  (4.19)  and  (4.20)),  of  n  variables  (the  poles  and  zeros).  This 
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problem  is  readily  solved  through  the  use  of  the  Fletcher- Powell  algorithm,  the 
details  of  which  can  be  found  in  Reference  20. 

A  program  that  performs  the  synthesis  of  recursive  digital  filters  using  the 
minimum  p-error  criterion  and  the  method  of  Fletcher  and  Powell  for  function 
minimization  is  contained  in  Reference  16.  Included  in  this  reference  are  an  exten- 
sive program  description,  input  requirements,  dimension  restrictions,  and  examples 
to  illustrate  how  the  program  is  used  for  various  types  of  filter  design.  Also  included 
is  a  copy  of  the  actual  program. 

The  reader  is  reminded,  however,  that  to  make  use  of  the  program,  the  transfer 
function  must  be  in  the  cascade  form  mentioned  previously.  The  program  output 
consists  of  an  argument  vector,  X,  and  a  gradient  vector,  G,  corresponding  to  the 
minimization  of  the  functions  Fl,  F2,  and  F3  which  are  used  by  the  program  for 
the  magnitude  approximation,  the  group  delay  approximation,  and  the  combined 
magnitude  and  group  delay  approximation,  respectively.  Based  on  X,  the  poles, 
zeros  and  coefficients  of  the  cascade  realization  of  the  filter  are  computed.  If  desired 
by  the  user,  the  frequency  response  is  also  given  [16]. 

D.    SUMMARY 

In  summary  then  it  is  appropriate  to  cite  a  recent  paper  by  Little  and  Gowdy 
[17],  that  evaluates  both  the  optimum  FIR  design  method  (that  employs  the  Remez 
exchange  algorithm),  and  the  minimum  p-error  method  used  in  IIR  design.  Their 
evaluation  specifically  investigates  convergence  problems  encountered  when  using 
these  iterative  techniques,  and  considers  the  design  of  lowpass,  highpass,  bandpass 
and  bandstop  filters. 
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The  following  is  a  brief  summary  of  the  more  important  points  mentioned  in 
this  article,  that  should  be  considered  when  using  these  iterative  techniques. 

1.  Minimum  p-Error  Design  Method 

a.  Advantages: 

•  Extremely  flexible  -  can  be  used  to  design  filters  with  arbitrary  magni- 
tude and/or  phase  characteristics,  as  opposed  to  being  limited  to  lowpass, 
highpass,  bandpass,  or  bandstop  designs. 

b.  Disadvantages: 

•  Nonconvergence  problems  due  to: 

o  a  poor  guess  for  the  initial  parameter  vector  X 
o  finite  wordlength 

•  Uses  large  amounts  of  CPU  (IBM  3081)  time  (up  to  2  minutes  for  higher 
order  filters,  greater  than  N  =  12) 

•  Requires  a  considerable  amount  of  input  to  specify  one  filter 

2.  Optimum  FIR  Filter  Design  Program 

a.  Advantages: 

•  The  program  is  easy  to  use 

•  Consumes  relatively  little  CPU  (IBM  3081)  time,  provided  the  filter  is 
not  of  extremely  high  order. 

b.  Disadvantages: 

•  Restricted  to  the  design  of  the  more  common  filter  types:  multi-band, 
bandpass,  Hilbert  transform,  and  differentiators 

•  Highpass  and  bandpass  filter  designs  were  poor  when  even  filter  lengths 
were  used,  but  good  for  odd  filter  lengths 

•  The  user  must  supply  an  FFT  program  to  obtain  frequency  response  data 
to  verify  the  filter  design.  Convergence  of  the  program  does  not  guarantee 
an  acceptable  design. 
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Keeping  these  considerations  in  mind,  both  programs  can  be  used  effectively 
to  design  a  wide  variety  of  FIR  and  IIR  filters.  It  is  suggested  that  if  a  designer 
wishes  to  use  either  of  these  programs  this  paper  is  a  valuable  reference  to  assist 
in  identifying  problems  that  may  be  encountered. 
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V.    CONCLUSION 

The  areas  of  recursive  (IIR)  and  nonrecursive  (FIR)  filter  design  have  been 
investigated,  while  the  more  predominant  design  methods  have  been  extensively- 
discussed  and  exemplified.  Additionally,  the  prevailing  computer-aided-design  al- 
gorithms (Remez  exchange  and  Fletcher-Powell)  were  also  presented. 

In  conclusion,  points  in  these  three  areas  that  merit  special  emphasis  are 
summarized  below: 

A.  RECURSIVE  FILTER  DESIGN 

•  The  desired  frequency  response  may  be  obtained  using  a  lower  order  filter 
than  if  a  nonrecursive  realization  were  used,  however,  filter  stability  must 
be  considered,  and  a  linear  phase  characteristic  is  not  guaranteed. 

•  The  traditional  design  methods  involving  Butterworth,  Chebyshev  or  el- 
liptic analog  prototypes,  and  the  bilinear  transformation  are  algebraically 
intensive.  The  direct  design  method  presented  eliminates  the  need  for 
determining  an  analog  prototype  and  for  prewarping,  thereby  reducing 
the  number  of  calculations  required. 

B.  NONRECURSIVE  FILTER  DESIGN 

•  Although  a  higher  order  filter  is  required  than  in  recursive  realizations  to 
obtain  the  same  frequency  response;  nonrecursive  filters  are  always  stable 
due  to  the  all  zero  nature  of  their  transfer  functions,  and  can  exhibit  a 
linear  phase  characteristic. 

•  The  filter  coefficients,  h(n),  may  be  determined  analytically  by  expanding 
the  desired  frequency  response  in  a  Fourier  series,  because  the  Fourier 
coefficients  correspond  to  the  filter  coefficients. 

•  Windows  may  be  used  to  eliminate  the  overshoot  in  the  frequency  re- 
sponse (Gibbs'  phenomenon)  caused  by  truncating  the  number  of  terms 
in  the  Fourier  series  to  a  finite  value.  However,  windowing  decreases  the 
sharpness  of  the  filter's  cutoff  region. 


165 


•  Frequency  sampling  can  be  used  when  an  analytical  expression  for  the 
desired  frequency  response  cannot  be  found.  The  IDFT  is  then  used  to 
determine  the  filter's  impulse  response  from  the  these  sample  values. 

C.    COMPUTER-AIDED  DESIGN 

•  For  FIR  filters,  programs  using  the  Remez  Exchange  algorithm  are  the 
most  popular,  while  for  IIR  filters,  the  Fletcher- Powell  algorithm  is  in 
predominate  use. 

•  CAD  is  especially  advantageous  in  the  design  of  extremely  high-order 
filters,  or  filters  with  arbitrary  frequency  response  characteristics. 

As  stated  in  the  introduction,  the  design  methods  presented  here  are  by  no 

means  the  only  ones  available  to  the  filter  design  engineer;  however,  d^e  ,o  the 

quantity  of  information  available,  one  is  easily  overwhelmed.   A  need  exists  for  a 

concise  source  of  information  on  the  popular  design  methods  available,  and  how 

they  are  used.  This  has  been  the  intent  of  this  thesis. 
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APPENDIX 
MAIN  PROGRAM 


c 

C  MAIN  PECGEA.1:  FIE  LINEAR  PHASE  FILTER  DESIGN  FECGEAM 

C 

C  AUTHORS:     JAMES    H.     HCCLELLAK 

C  DEPAETSSHT     Of     ELECTRICAL    ENGINEERING    AND    CCMPUTEF    SCIENCE 

C  MASSACHUSETTS    IBSlITlllI    CF    TECENCLCGY 

C  CAMBRIDGE,     MASS.     0^1  3£ 

C 

C  THOMAS  H.  PARKS 

C  DEPARTMENT  OF  ELECTRICAL  ENGINEEEING 

C  RILE  UNIVERSITY 

C  HOUSTON,  TEXAS  77001 

C 

C  LAWRENCE  E-  RABINEE 

C  BELL  LABORATORIES 

C  MORBAY  HILL,  NEK  JERSEY  07S74 

C 

C  INPDT: 

C  NFILT--  FILTER  LENGTH 

C  JTYPE —  TYPE  CF  FILTER 

C  1  =  MULTIPLE  PASSbAND/STCPBAKE  PILIEE 

C  2  =  DIJFEEENTIATCF 

C  3  =  HUBERT  TRANSFORM  EILTE2 

C  NBANDS —  NUMBER  OJ  EANDS 

C  LZL1D —  GRID  DENSITY,  illl  EE  SET  TC  16  UNLESS 

C  SPECIFIED  OTHERWISE  BY  A  POSI1IVE  CCNSTANI. 

C 

C  EDGE(2»NBANDS)  —  3ANDEDGE  AERAY,  LOWER  AND  UFPEE  EDGES  FCF.  EACH  B^ 

C  i'llb  h    MAXIfiCC  Of  10  EANDS. 

C   FX (NBANDS) —  DESIRED  FUNCTION  AEEAY  (CE  DESIRED  SLCEE  IF  A 
C  DIFFERENTIATOR)  tCE  EACH  BAND. 

C   WTX(NoANDS) —  WEIGHT  FUNCTION  ARRAY  IN  EACH  BAND.   FOE  A 

C  DIFFERENTIATOR,  TEE  WEIGHT  FUHCTICii  IS  INVERSELY 

C  PROPORTIONAL  TC  I. 

C 

C   SAMPLE  INPUT  DATA  SETUP: 

C  32,1,3,0 

C  0.0,0.1,0.2,0.35 

C  0.-25,0.5 

C  0.0,1.0.0.0 

C  10.u,1.0,1u.0 

C  THIS  DATA  SPECIFIES  A  LENGTH  32  BANDPASS  FILTER  WITH 

C  STCPBANDS  0  TC  0.1  AND  C.-25  TC  0.5,  AND  PASSBAND  FBCfl 

C  0.<<  10  0.35  WITH  WEIGHTING  CF  10  IN  THE  STOPBAuDS  AND  1 

C  IN  THE  PASSBAND.   THE  GPID  DENSITY  DEFAULTS  TO  1o. 

C  THIS  IS  THE  FILTER  IN  FIGURE  10. 

C 

C  THE  FOLLOfllNG  INPUT  DATA  SPECIFIES  A  LENGTH  32  FULLEAND 

C  DIFFERENTIATOR  KITH  SLOPE  1  AND  WEIGHTING  CF  1/F. 

C  THE  GRID  DENSITY  WILL  EE  SET  TC  20. 

C  32.2,1,20 

C  0,0.5 

C  1.0 

C  1.0 

C 

c 

COMMON  PI2,AD,DEV,X, Y, GRID ,DES, WT, ALPHA, I EXT , NFC NS, NG RID 
COMfiCN  /OOPS/MTEr,ICUT 

DIMENSION  I  EXT  (bo)  ,AD(66)  ,  ALPHA  (6  6)  ,  X  (b6)  ,  Y  (66) 
DIMENSION  H  (bo) 

DIMENSION  DES (10U5)  .GRID (104  5)  , WT (104  5) 
DIMENSION  zlGil  (20)  ,FX  ( 1  0)  ,  KIX  ( 10)  ,DEVIAI  (10) 
DOOdLE  PEECISION  FI2.PI 
DOUBLE  PEECISION  AD,D£V,X,Y 
DOUBLE  PRECISION  GEE,D 
INItGEE  BD1,fcD2,EE3,BD« 

DAI  A  BD  1,EE^,  ED  J,  BD-+/1HB,  1flA,  1BN,  1HD/ 
C      lNFUI  =  I1iiACH  (1) 
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C  ICDI=I1BACH(2) 

P1=4.0*DAIAK  (T.ODO) 

FI2=2.QDOO*P1 
C 

C       THE    PECGEAB    IS    SET    OP    FOE    A    BAXIBUB    LffiGTE    CF    128,     EUT 
C       THIS    UPPE*    Llilll    CAN    BE    CHANGED    EY    EEEIBENSICNING    THE 
C      AREAYS    1EXX.    AC,    ALPHA-     X,     X,    H    TO    BE    NFBAI/2    ♦    2. 
C       THE    ABtvAYS    DES,     GBID,     AND    KT    BUST    DIEENSICNED 
C       1t>  (NFBAX/2    ♦    2)  . 

NFBAX=123 
100  CONTINUE 
JIYPI=0 
C 

C   PHCGEAB  INPUT  SECIICN 
C 

BEAD  (4,110)  NFI1T,JTYPE,NEANES,LGEID 
If  (NFILT.EQ.O)  S1CP 
110  FGRBAT  (4  (Id,3X) ) 

IF  (6iFIlI.LE.NFaAX.ANL.NFIIT.GI. 3)  GO  TC  115 
CALL  EEECE 
STCP 
115  IF (NBANJS.LZ.O)  NcANDS=1 

C   GEID  DENSITY  IS  ASS02ED  TO  BE  16  UNLESS  SFECIFIEC 
C   CIHZEhlSE 


C 


IF  (LGEID.LE.O)  IGEID=1o 
JB=2*NBANDS 

:AD(4,1L0)   (EDGE  (J)  ,J=1,JE) 

[FX  (J)  ,J=1,NBANDS) 


4,120) 
BEAD  (4,  120)   (WIX  JJ)  ,J=1,  lib  AND  3) 
FCB3A1  {hI15.4) 
IF  (JIYPE.G1.0.ANE.JTYPE.LE.3)  GO  TO  125 
CALL  EEhCB 
STOP 
125  NEG-1 

IF  (JTYFL.EC..1)  NEG=0 

NCLD=NFILI/2 

«UDD=NFILT-2*NCDD 

NFCNS=N5ILI/2 

If  (NCDE.EQ. 1.ANE.HEG.EC-0)  NFCNS  =  Nf CNS  + 1 

C   SET  UP  THE  DENSE  GEID.   THE  KDHcEfi  Of  POINTS  IN  IKE  GBID 
C   IS  (FILTER  LENGTH  +  1)*GEIE  EENSITY/2 


GBID  (1)=ZDGE  (1) 

DELf=LG5IE»NFCNS 
DELf=0. 5/DELf 


F  (NEG.EC.O)  GO  TC 
^ (1) -LI.Dr 
135  CONIINOE 


IF (EDGE  (T) -LI.DELF)  oRID(1)=DELF 


G=1 

I  =  1 

IBANL-1 
140  FUF=EDGE(1+1) 
145  IEBP=G£ID(J) 

C   CALCULATE  THE  DE3IEEE  BiGNITUDE  EESPCNSE  ANE  THE  WEIGHT 

C   FUNCTION  CN  IHE  GRIE 

C 

DES  (J)=EFF  (TEBP,F X , XTX , LB A ND ,JTYPE) 

KT  (J)=WA1E  <IEBP,FX,i.IX,LaAND#JTYPEJ 

GEID  (J) =TEBP+DELF 

II  (GEID  (J) -GI. FOP)  GC  TC  150 
GC  TC  145 

150    GEID  (J-1) =FUP 

DES  (J-1)=EFF  (FUP,FX,WTX, LEAKE, JTYPE) 
WT  (J-1)=iATE  <iUF,FX,WTX,LfcANE,JTYPE 
LBAND  =  LBAND«-1 
L  =  I  +  2 
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IF  (LHAND.GT.NEANES)  GO  1C  160 

GEIC  (J)  =  EEG£  (1) 

GC  IC  -WO 
160  NGEID=J-1 

IF  (NtG.NE. NCED)  GC  TC  165 

IF  (SLID  (NGEID)  .GT.  (0.5-EELF)  )  NGEID=KGEID- 1 
165  C0M1IHDE 
C 

C   SET  DP  A  NEU  APP5CXIMATICN  PBCELEK  WHICH  IS    EQUIVALENT 
C   TC  1HE  GBIGINAl  FECEIEH 
C 


IF  (NEG)  170.170,180 
170  IF  (NCED.EC-  1)  GC  IC  200 

EC  175  JM.NGHIE 

CRANGE=DCCS(?I*GEICJJ) ) 

DES  (J)  =  i)ES  (J) /CHANGE 
175  WT  (J)=K1  JJ)  *CriANGE 


GO  Tl 
IF  (N( 
DC  145  J=T,NGEIS 


180  IF  (NCDD.EC.1)  GO  TC  190 


CHANGS=DSIN  (F1*GEID  (  ")  ) 

DES  (J) =DES  (J) /CHANGE 
185  WT  (J) =WTJJ) "CHANGE 

GC  IC  200 
190  DO  195  J=1, NGEID 

CEANGE=DSIN  (?I2*G5ID  (J)  ) 

EES (J) =DES   (J) /CHANGE 
195  »'l  (J)  =UT  (J) 'CHANGE 

C       INITIAL    GUESS    FOE    IliE    EXTEZflAL    FEEQUENCIES — ECUALLJ 
C       SPACED    ALONG    THE    GFII 


C 


200  TEHr=FLOAT  (NGEIE-1) /FLOAT  (NFCNS) 

DO  210  J=1, NFCNS 

XI=J-1 
210  IEXI  (J)=XT*IEHP*1.0 

IZXI  (KFCiiS+1)  =NGEIE 

NH1=NFCNS-1 

NZ=NFCNS*1 
C 

C   CAL1  THE  &EJ1EZ  EXCHANGE  ALGORITHM  TO  EC  IHZ  APPECX IHAI1CS 
C   PEGEIEM 
C 

CALL  SZKEZ 
C 
C   CALCULATE  THE  IflPULSE  KZSPONSZ. 


C 


IF(.NEG)  3C0, 300,320 
300  IF  (NODE. E¥. 5)  GC  IC  310 

EC  jC5  J=1,NH1 

NZ5J=N2-J 
305  h  (J) =0.5*ALPHA  (NZflJ) 

E(NFCNS^=ALPHA  1) 

Gd  IC  350 
310  H<1) =0.25»ALPHA (NFCNS) 

EO  315  J=2, NB1 

N2HJ=NZ-J 

NF2J=NFCNS+2-J 
315  H(J) =0.25»  (ALPHA  (HZMJ)  ♦ALPHA  (NF2J)  ) 

h  (NFCNS)  =0.5» ALPHA  (1)  «-0.25» ALPHA  (2 

GC  IC  350 
320  IF  (NCDD.EC-0)  GC  IC  330 

H (1) =G.25*ALFHA  (NFCNS) 

H(2[  =  0.25*ALFHA  (NH1) 

EC  325  J=3,NM1 

NZHJ=NZ-J 

NFJJ=NFCNS+3-J 
325  H  (J) =0.25*  (ALPHA  (HZttJ) -ALPHA (NF3J) ) 

H  <NFCNS)=0.5»ALPh*  ( 1) -0. 25»ALFHA  (3 

h  (NZ)=0.0 

GC  IC  350 
330  H  (1) =0.25*ALFHA (NFCNS) 
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DC    335    J=2,NM1 
NZflJ=KZ-J 
NF2J=NfCNS+2-J 
335    H  (J)  =0.25*  (ALPHA  (NZflJ)  -ALPHA  (NF2J) ) 
fa  (NFCNS)  =0.5*  ALP  HA  ( 1j -0.  25*  ALPHA  (2) 

C       PHOGBAfl    OUTPUT    SECIION. 

C 


350  'JEIIE  (8,360) 
360  FCEflAT(1H1,  70  (1 H*)  //  15X,  29HFIKITE  IMPULSE  EESPCNSE  (FIB)/ 
1 13X,34HLINZAE  PHASE  DIGITAL  F1LTIE  DESIGN/ 
217X,24HEEBEZ  EXCHANGE  ALGCEITHB/) 

IF  (JTI?£.E;.1)  KBIIE(8,3o5) 
365  FOEflAI  (221,  15HBAN2PASS  i'lLIEF./) 

If  (JIYPE.EC..2)  WEIIE]d,370) 
370  POEHA1  (221, 14HDIFFEEENIIAICE/) 

IF  (JTYPE.Ei.  J)  iEIIE  (8,375) 
375  F0E3AT  (20X, 19HHILBEEI  TEANSF CEHEB/) 

WEIIEj6\j78)  NFILI 
378  FOEMAI (20X, 16HFILIEE  LENGTH  =  ,13/) 

VBIIEJ6, j80) 

380  POBflAI  (151,288*****  IflPOL.SE  RESPONSE  ****») 
DO  381  J=1, NFCNS 

K=NFILT*1-J 

IF  (NZG.Ew.O)  WBII£(b,382)  u.B(J),F. 

If  (KEG. Ew-1)  «HIIE{a,383)  J,H(Jj,K 

381  COMIIIDE 

382  FOEflAT  (13X,2HH  (,I2,4H)  =  ,E15.3,5H  =  H(,I3,1U)) 
38J  FORMAT  (13X,2HH (,I2.4H)  =  ,E15.8,6H  =  -a(,I3,1H)) 

IF  (NEG.Ei.1.  AND.JiCiJ.ES.  1)  SSIIZ(b,334)  NI 

384  FCRBAI  (13X, 2HH (, 12, 8H)  =   CO) 
CO  450  K=1,NBANDS,4 

KDF=K+3 

IF  (KUP.GI.NBANDS)  KU?=N3ANDS 

KBIT E (8,365)   (dE 1 , ED2.ED3, ED4, J ,J=K, KUP) 

385  FOEBAT  </24X,4  (4A1,Ij,7X)  ) 

utilli  (o,390)   [ECGt  (2*J-li,J=K,K0P) 
390  FORMAT  12X,  15HLC'.EE  BAND  EDGE, 5f 14. 7) 

•  SITE  (0,395).  (EEGc  (2»J)  ,J=K,K0P) 
395  FCKnal  (2X,l3Hu?PEF  BANE  EDGE, 5F1U. 7) 

IF  (JT1PE.NE.2)  5EITE(8.4J0)   (f X  [ J)  , J=K , KUP) 
400  FCEBAI  (2X,13HUES1RED  V ALUE ,2X, 5* 1 4 . 7) 

IF  (JTYPE.EC.2)  WEITE  (8.4J5)   (F X  ( J)  , J  =  K, KUP) 
405  FOEBAT  (2X, 1 3HDESIEED  SLOPE ,2X, 5F 14. 7) 

WBIIE(e,4lO)   (VIXJJ) ,J=K,KUP) 
410  FOBflAT  (2X,9EV£IGEIING,6X,5F14.7) 

CO  420  J=K,KU£ 
420  EEVIAT  (J)=D£V/KIX(J) 

WRITE (6,425)   (DEVI AT  (J)  ,J=K-KOF) 
425  FCEBAT  (2X , SBDEV1AIIOK, o£,bF 14. 7) 

IF  (JIYPE.NE.  1)  GC  TO  450 

EC  430  J=K,KUP 
430  DEVIAT  (J) =20.0*ALCG 1 0 (DEVI AT (J)  *FX  (J) ) 

fcEIIE  (6,435)   (DEVIAT  (J)  ,J  =  fc,KUP) 
435  FCEKAI  (2X, 15HDEVIA1I0N  IN  EE,5F14.7) 
450  CONTINUE 

DO    45^:    J=1,N2 

IX=IEXT  (J) 
452    G£ID(J)=GEID(IX) 

5.EITE  (6,^55)      (GRID  (J)  ,J=1,NZ1 
455    FOBBAT  (/2X.h7HEXT££KAL    FEEwUENCIES  — H  AXIMA    OF    THE    EEECE    CDEVE/ 
1     (21,5*12.7)) 

WRITE  (6,460[ 
4o0    FOE  BAT  f/1X,70  (1H*J/1H1) 

GC  TC  100 

END 
C 

C 

C  FUNCTION:  EFF 

C    FUNCTIOr.  TO  CALCULATE  THE  DESIRED  KAGNIIUDE  EESPOHSE 

C    AS  A  FUNCTION  Of  FREQUENCY. 

C    AN  AEBITEARY  FUNCTION  OF  FREQUENCY  CAN  BE 
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C  APPROXIMATED  IF  THE  USER  REPLACES  THIS  FUNCTION 

C  WI1H  IHE  APPROPRIATE  CODE  10  EVALUATE  THE  IDEAL 

C  MAGNITUDE.   N01E  THAI  TflE  PARAMETER  FEEQ  IS  THE 

C  VALUE  Of  NORMALIZED  FREQUENCY  NEEDEE  FCR  EVALUATION. 

C 

FUNCIION  EFF  (FREC  ,  FX,  WTI  ,LEA  1JD,  JTYPE) 

ElEENSICN  FX  lb)   ,WIX  (5) 

If  (JTYFE.EQ.2)  GC  TO  1 

EFF=FX  (LBANE) 

RETURN 
1  EFr=FX (LBANE) *FEEQ 

RETURN 

ENE 
C 

C 

C  FUNCTION:  MATE 

C    FUNCTION  TO  CALCULATE  THE  WEIGHT  FUNCTION  AS  A  FDNCIICN 

C    Of  FREQUENCY.   SIMILAR  IC  IHE  FUNCTION  EFF,  THIS  FUNCTION  CAN 

C    BE  EEFLACEE  til  A  DSER-WRITIEN  ROUTINE  TC  CALCULATE  ANX 

C    DESIRED  BSIGHIIKG  FUNCTION. 

C 

c 

FUNCTION    KATE  (FEEQ , FX, WTX, LB AN £, JT XPE) 

EIMENSION    FXjb)  ,»IX  (5) 

IE  (JIXFE.EC.2)  GC  TO  1 

UAIE=K1X(LEAND) 

RETURN 

1  IF  (FX  (LBAND) .LI.0.C001)     GC    TC    2 
UAIE=WIX  (LBANDJ/EEEQ 

RETURN 

2  WAIE  =  WIX  (LBANE) 
RETURN 

ENE 
C 

c 

C  SDBECUIINE;  EBECF 

C    THIS  FCUTINE  WHITES  AN  ERROB  MESSAGZ  13  AN 

C    EEBOfi  HAS  BEEN  EEIECTEB  IN  TEE  INPUT  DATA. 

C 

SUBROUTINE  ERROE 

coaaoii  /ccps/niizr,iohi 

■  RUE  (8,1) 
1     FORKAT  (44H    ««»»*»*»»»**    ERBCS    IN    INPUT    CAIA    «»«»»»**»•) 

RETURN 
ENE 
C 

C  SU3BCCTINE:  RZMEZ 

C    THIS  SUBROUTINE  IMPLEMENTS  IHE  RZMEZ  EXCHANGE  ALGORITHM 

C    FOR  THE  WEIGHTED  CHZtiYSEEV  AEPECXIM AT 1C N  OF  A  CONTINUOUS 

C    FUNCTION  -1Tb    A  SUA  OF  CCSINES.   INPUTS  TC  THE  SUBROUTINE 

C    ARE  A  DENSE  GBIE  WHICH  REPLACES  THE  FREQUENCY  AXIS,  THE 

C    DESIRED  JUNCTION  ON  THIS  GRID,  THE  WEIGHT  FUNCTION  CK  THE 

C    GRID,  THE  NUMBER  CF  COSINES,  AND  AN  INITIAL  GUESS  OF  IHE 

C    EXTREMAL  FREQUENCIES.   IHE  PROGRAM  MINIMIZES  THE  CHEcISBEV 

C    ERROR  EY  DETERMINING  THE  EEST  LOCATION  CF  THE  EXTREMAL 

C    FREQUENCIES  (POINTS  OF  MAXIMUM  ERFOB)  AND  1HSN  CALCULATES 

C    THE  COEFFICIENTS  CF  THE  BES1  APPROX IMATICN. 

C— — — — ———————— — — — — — — — — —  _ — _ _ _— __— __ — — — _ 

c 

SUBROUTINE  EEMEZ 

COMMON  PI2,AB,DEV,X, 2. GRID ,BES , KT , ALP  HA, IEXT, NFCNS, NGBID 

COMMON  /COPS/NITER, ICUI 

EI MENS ION  IEXT(oo[,AB(66)  ,ALPhA(6o)  ,  X  JbG)  ,  Y  (06) 

DIMENSION  DES  (104  3)  ,GEI£(1045)  ,  Ml  (1045) 

DIMENSION  A  (oo)  ,F  (to)  , Q  (o5) 

DOUBLE  PRECISION  t 12, DNDM, BEEN , DTEMP, A, E ,Q 

ECU3LE  PRECISION  ER,EAK 

DUDBLE  PRECISION  AD,DEV,X,X 

DOUBLE  PRECISION  GEE,D 
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C   THE  PSOGEAd  ALLOWS  A    8AXIHUH  ifOHEEB  CF  ISIBAIIOKS  OF  25 

ITEflAX=25 

EEVL=-1.0 

N2=KFCUS+1 

KZZ=NfCHS+2 

KIIEf=0 
100  COKIIHDE 

IEXI  <NZZ)  =KGEID*1 

NITEF=NI1EE*1 

If  (NITER. GT.IIE8AX)  GC  TC  40C 

EC  110  J=1,NZ 

JXT=IEXI <J) 

ETEHr=GEIE (JXT) 

DT£flE=ECCS (EI£flP*PI2) 
110  XjJ)=LIEHP 

JEI  = JSPC«S-1)/15+1 

DO  120  J=1,NZ 
120  ADJJ)=D  <J,NZ,JET) 

EK0B=0-0 

EE£N=0.0 

K=1 

EC  130  J=1,MZ 

L  =  IcXl  <J) 

EI£KP=Afl (J)*D£S (L) 

DNuH  =  DNU£+DIE.lF 

Dl£flf  =  FLOAI  <K)  *AE  (JJ/KT  (L) 

E£EN=E£EN*EI£flP 

130  K=-K 
E£V=DKU£/DDEfl 
kEITE <8, 131)  BEV 

131  F0E3AI <1X,12HD£VIAII0K  =  ,F12.9) 
NU=1 

If  (EEV.GT.0.0)     N0  =  -1 

EEV=-f~CAI  <KU    »E£V 

K  =  K0 

EO     14C    J=1,NZ 

L=IIZ1 <J) 

El£."-f  =  riOAT  (K)  *EEV/il  (I) 

i'  i0)  =DES  (I)  *£IZfir 
140  K  =  -K 

IF-iDEV.GT.DEVi)  GC  TC  150 

CALL  OUCE 

GC  TC  400 
150  EEV1=D£V 

JCHNG£=0 

E1=IIXT{1) 

KN2  =  I£XT  (NZ) 

KLCW=0 

N0T=-NU 

J=1 
C 

C   SEAECE  FOE  THE  EXTEEflAL  FEEQUFNCIES  CF  IHE  BEST 
C   A?P£CXIHAIICN 
C 

200  IF  (J.EC..HZZ)  YNZ=CCHP 

If  (J.GE. NZ2[  GO  TC  300 

KUP^IEXT (J*1) 

L  =  IIX1  (J)  +1 

KUI=-NU1 

If<J.£C-2)  J1=CCHP 

COCP=DEV 

IF  JL.GE  .F.0F1  GC  TC  220 

EEE  =  G££  (L.NZ) 

£EE=  (EEE-EEJ  (L)  ) *HI  (L) 

ET£tr=FLCAl  (NUT) *EBE-CCHP 

IF  (ETEMP.LE.G.O   GC  1C  220 

COBP=FLGAT  (NUT) *EBB 
210  1*1*1 

IF  (L.GE.KUP1     GO    TC    215 

EBJ=GEE  (L,NZ) 
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E5E=  (EEF-DES (l) ) *II  11) 

EI£»i£  =  FLOAl  (NUT  »EE*-CC3P 

IF  (EIEBF.LE.O.O)  GC  TC  215 

COhP  =  FLOkl   (NUI) *EEfl 

GC  IC  210 
215  IEXT(J)=1-1 

J  =  J+1 

KICli  =  I-1 

JCHNGZ  =  JCHNGF,+  1 

GC  TC  200 
220  1=1-1 
22D  1=1-1 

IF  (1.1E.KLCW)  GC  1C  250 

EEF=GEE (1,NZ) 

£EF= (ZRF-DES  (I)  )  *WI  (L) 

ETZMP=FL0A1  (NOT)  »E£*-CCMP 

IE  (DlZaF.GT.O.O)  GC  1C  230 

IF  (JCfiNGE.IE.O)  GC  TC  225 

GC  IC  2b0 
230  C0f1P=F10AI  (NUI)  *£EE 
235  1=1-1 

If  <1.1£. K1CW)  GC  TC  240 

EEE=GEE  (I,K2) 

ERB=  (Eah-DES (L) ) *WT  (1) 

EIZHP=FLCA1  (NUI)  *EE?-CCMP 

IF  (ETEKP.LZ.G.O)  GC  TC  240 

CUaP=FLCAI  (NUI)  »EEK 

GC  TC  235 
240  KlCk=lZXT(J) 

IZX1  (J) =  L+1 

J  =  J*1 

JCHNGF=JCHKGF+1 

GC    IC    2v>0 
250    1=IZZI (J) +1 

IE  (JCH»«E.GI.O)     GC    TC    215 
255    1=1-1 

IF  (L.GE.K0P1     GC    TC    2oO 

£EE=GZE  (L,KZ) 

Lkh=  (EEE-DES(L) ) »KI (L) 

ETZnP=FLCAI  (2.UI    *EEE-CCf1P 

IF jEIEflr.li.C.O)     GC    10    255 

CCKP=FLCAI  (NUI) *£EE 

GC    TC    210 
260    KlCi.=IZXl  (J) 

J  =  J+1 

GC  IC  200 
300  IF  (J.GI.HZZ)  GC  IC  320 

IF  (K1.GI.IIXT  (1) )  K1  =  IZXT(1) 

1FJKNZ.LI.IEXT  (NZ) )  KNZ=IZ2T  (NZ) 

N0T1=N0T 

KUI  =  -i<U 

1=0 

KUP=K1 

COHP=YNZ* (1.00001) 

1UCK=1 
310  1=1*1 

IF  (l.GE.EUF)  GC  TC  315 

EKF=G£E  (L,NZ) 

EEE  =  (lafi-DES(L)  )  *WT  (L) 

ETEBF=FLCAT  (NU I) *E RE-CCflP 

IF  (DIZHP.LE.Q.O   GC  IC  310 

CCnF=FLOAl  (NUT) *£EE 

J=N2Z 

GC  TC  210 
315  LUCK=6 

GC  TC  325 
320  IF  (LUCK. Gl. 9)  GC  TO  350 

IF  (CCKE.G2.il)  Y1=C0HP 

K1=IEiT  (NZ2) 
325  I=NGEID+1 

KLC»=KNZ 

NU1=-NUT1 
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C0HP=I1»  (1.00J01) 
330  L=L-1 

IF  (L. IE. KICK)  GC  IC  j40 

EEr=GEr  (L, HZ) 

£EB={ZSi-EES  [111  *HT  (L) 

CTEfl£=FLCAl  (NUI)  *£EE-CC!1P 

IF <D1EHP.LE.0.0)     GC    IC    330 

J*K2Z 

COfll*PIOAI  (NUI) *E5B 

lOCK-IDCK+10 

GC    IC    235 
340    IF  (LUCK-EC. 6}     GC    TO    370 

DC    345    J=1,HFCH2 

NZ2flJ=NZ2-J 

K2tiJ  =  NZ-J 
345    IEXT  (NZ2KJ) =IEXT  <NZMJ) 

1EXIJ1)=K1 

GO  IC  TOO 
350  KN  =  IEXT  (N2Z) 

DO  360  J=1,S1CMS 
360  IEXIlJ)=IEII  (J  +  1) 

IEXI  (KZ)  =  KN 

GC  IC  100 
370  IF  (JCflnGE.GT.O)  GC  IC  100 

C   CALCULAIICN  OF  THE  COEFFICIENTS  OF  THE  cE3I  APPECXIM AIIC N 

C   UJlJiG  I  HI  INVtSSE  EISCEZTE  FCOEIEB  TBAKSFCLfl 

C 

400  COBIIKOE 

Nfll=NFCNS-1 

FSH=1. JE-Ofc 

GIEMt=GEICJ1) 

XIN22)=-i.O 

C»a2*NFCMS-1 

EELF=1.0/CN 

L=1_ 

IF  (GFID  {1) .LT. 0.01. AND. GEIE(NGBID)  .SI.C.49)     KKK  =  1 
IF  (NFCNS.LE.3)     KKK=1 
IF  (KKK.EQ.  1)     GC    IC    <*05 
DIEfiF»DCUSfPI2*G5IDMn 
EKUflsJJCCS  (5I2*GiIE  (wGEID)  ) 
AA=2.0/  fDIEflP-DKUS) 
EB  =  - (CIEH?*DNOaj/  (CIEBP-DN Ofl) 
405    COKllNUE 

DC    43C    J=1,NFCMS 

FT=J-1 

FI=F1*DELF 

XT=ECCS  (PI2»FI) 

IP  (KKK.2S.1)     GO    IC    410 

11= <X1-B£)/AA 

XI1  =  SgBT  (1.0-XT*XI) 

FT=A1AK2{IT1#XT)/EI2 


410    XI=X(L) 

IF  (II.GT.XE1     Gl 

IF  (  <XE-XI) .Il.FSb)     GO    10    4  15 


IF  (II.GT.XE1     GC    IC    420 


1  =  1+1 

GC    2C    410 
415    A(J)=X{L) 

GO    IC    425 
420    IF (jXT-XE) -II. FSH)     GO    IC    415 

GEID<1)=FI 

A(J)=GEE(1.NZ) 
425    CONT1N0E 

IF  (L.GI.  1)     I  =  L-1 
430    CONIIKOE 

GEID  (1)=GIZKP 

ED£N=PI2/CN 

CO    510    J=1,NFCNS 

E1EMF=0.0 

ENUfl=J-1 

EhUfl=DNUfl»EDEN 
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If  (Nfl1.LT. 1)  30  IC  505 

EC  500  K=1,N31 

EAK  =  A  (K+1) 

EK  =  K 
500  DI£KF=DTEHF+DAK*ECCS  (DNUfl-EK) 
505  DlE«£«2.0*£ISflP*A(l) 
510  ALPHA JJ) =EIErtP 
-    550 


EC  550  J=2,NFCNS 

ALPHA  (J) =2.0*ALFHA (J)  /CK 

ALPHA  (1) =ALPHA (1) /CN 

IP  |KKK.£C.  1)  30  1C  545 

?  (1)=2.0*ALFHA  (til  CNS)  *Bri+AlPEA  (Nfll) 

?  <2) =2.0-_A*ALFHA jNFCNS) 


•)  =2.0*AA*ALFHA  (NFCNS) 
)=  ALPHA  (MFCHS-i) -ALPHA  (JiFCSi 

540  J=2,Nfl1 


If  (J.LI. Nfll)  GO  IC  515 

AA=0o5*AA 

£B=0.5»BB 
515  CONTINUE 

P  <J+1)=0.0 

EC  520  K=1,J 

A{K)  =P  (K) 
520  P  (K) =2-0»E£*A (K) 

?  (2)  =P  12)  «-A  (1)  »2.0*AA 

JH1-J-1 

EC  525  K=1,Jfl1 
525  FJK)»P(K)+C.(K)+AA*A(K+1) 

JP  1=J  + 1 

EC  530  K=3,JP1 
530  PJK)  =P  (K)  «-AA*A  (K-1) 

If  (J.2C.  NS1)  SO  TC  540 

EC  335  K=1,J 
535  ;  (K)  =-A  (K) 

NF1J=Nf_NS-1-J 

Q (1)  =C (1) *ALFHA  (Nf 1J) 
540  CONTINUE 

DC  54j  J=1,NfCNS 
343  ALPHA  (J) =P  <J) 
545  CONTINUE 

IF (NFCNS.GT.3)  F.EIUBN 

ALPH*  (NFCNS+1)=C.0 

ALPHA  NFCN5+2  =0.0 

HEIUEN 

END 


C  FUNCTION:  E 

C    FUNCTION  TC  CALCULATE  THE  LAGF.ANGE  I NTEEFCLATICN 

C    COEFFICIENTS  FOB  USE  In  IHE  FUNCTION  GIL. 

C— ——--——— — — — — ________ — — - — ______________________ 

c 

DOUBLE  PRECISION  FUNCTION  E(K,N,fl) 

COB HON  P12,AE,DEV,X, i  ,  SRIE  ,DE3,  WT  ,  ALFliA,  IE  XT  ,  NF 

EISENSIGN  IEXT  (bo)  ,A£  (bo J  .  ALPHA  (*o)  ,  X  (bo)  ,Z  (b6) 

EIHENSICN  DES  (104b)  , GEID ( 1 04  5)  ,  31  ( 104  5) 

EOUBLE  PRECISION  A£,DEV,X,X 

BOUELE  PRECISION  C 

EOOfilE  PEECISION  PI2 

E=1.0 

C  =  XJK) 

EC  j  1=1, fl 

EO  2  J=L,N,fl 

IFjJ-K)  1,2,1 

1  D  =  i-C*b-(w-X(J)) 

2  CONTINUE 

3  CONIINUE 
D=1.0/B 
EETUEN 
ENE 

C 

C  FUNCIICN;  GEE 
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C    FUNCTION  TC  EVALUATE  TBE  FEigOF.NCT  RESPONSE  USING  IHE 

C    LAGLAUGE  INTEEPCLATICN  FCBKULA  IK  THE  EAEYCENIEIC  FCSH 

C ■ 

C 

DOOBXE  PRECISION  FUNCTION  GIE(K.N) 

COHHCfl  PI2,AE,DIV,I,I,GEIE,DES,»If ALEH A,IEXI , NFCNS, NGEID 

EI3ENSI0N  IEXI  (60)  ,AL  (bo)  f  ALEHa(bC)  .X  Jot)  ,Y  (00 J 

r~a£KSICN  VlSiiQuh) ,SEID(1G45) , 81 <1045) 

CUJBLE  PRECISION  E,C,D,XF 

DCUHLE  PRECISION  112 

DOUBLE  PHECISION  A£,DEV,X,Y 

E  =  C.O 

XF=GBID  (K) 

XF=DOCS  EI2*XF) 

E=C.O 

DC  1  J=1,N 

C  =  XF-X(J) 

C  =  AD(J)/C 

E=D*C 
1  E  =  E*C*Y(J) 

GEE=E/D 

SEIDEK 

END 
C 

C  SCBECUIIKI:  CUCH 

C  *?.ITES  AK  EEROF.  MESSAGE  WHIN  THE  ALGCRITHE  FAILS  TC 

C  CCNVEEGE.   1HERE  SEEM  TO  BE  2.C  CONEIIICMS  UN0E5  WHICH  • 

C  THE  ALGCEITHtl  FAILS  10  CCNVEEGE:   (1)  THE  INITIAL 

C  GUESS  FCE  THE  EXIEE3AL  FEECUZNC1ES  IS  SC  PCCE  THAT 

C  THE  EXCHANGE  ITERATION  CANNC1  GET  SIAETEE,  CF. 

C  12)     NZAE  THE  IE5HINAIICN  CF  A  CCIifiECI  EESIGK, 

C  THE  DEVIATION  DECREASES  DUE  IC  BOUNDING  EIECF.S 

C  A;JD  THE  ESCGEA"  SICFS.   IN  THIS  LA1TEB  CASE  THE 

C  FILIEZ  DESIGN  IS  FRCDAHLY  ACCEEIASLE,  EUT  SHCOLD 

C  BE  CHECKED  BY"  COBEUTING  A  FEEC.OZNCY  RESECKSZ- 

C 

SUBBOUIINE  CUCH 

CCH«CN  /CCES/NITIE,ICDI 

SiRIIE  {6,  1}  NITEE 
1  F05HA.  i-^L  »»»»»»»**»*»  FAILUEE  TC  CCNVEEGE  *«»»»*»***/ 
1**1H0?EOliAELZ  CAUSE  IS  flACBINE  HCUKDING  EEEOE/ 
223H0NUUBEE  CF  IIEEATICNS  =,14/ 
3J9H0IF  THE  NUH3EE  OF  ITERATIONS  EXCEEDS  i,/ 
4o2H0THE  EZSIGN  SAX  BE  CORRECT,  BUT  SHOULD  BZ  VEEIFIEE  CIIH  Afc  FF: 

EEIOEN 

END 
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